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Abstract 


Gel  propulsion  systems  combine  the  best  characteristics  of  solid  and  liquid  propellants. 
The  gel  system  stores  like  solid  propellant,  but  flows  like  a  liquid  when  pressurized, 
enabling  throttle  and  restart  capability  similar  to  liquid  propellants.  In  addition,  gels  have 
a  lower  vapor  pressure  compared  to  their  liquid  counterparts.  Consequently,  the  storage 
containers  holding  them  can  be  much  lighter. 

An  enabling  technology  for  the  further  advancement  of  gelled  propulsion  is  the 
development  of  tools  to  render  the  synthesis  of  propellant  formulations  more  systematic. 
In  Phase  I,  various  models  were  developed  to  support  the  project  objectives.  Given  the 
multi-scale  scope  of  the  problem,  first  quantum  mechanical  and  Grand  Canonical  Monte 
Carlo  (GCMC)  simulations  were  used.  The  results  of  these  simulations  were  then  fed 
into  Stokesian  models  to  compute  stress  characteristics.  Finally,  Brownian  and 
continuum  models  were  used  for  syneresis  studies. 

In  Phase  II,  Quantitative  Structure-Property  Relationship  (QSPR)  techniques,  in  the  form 
of  a  feed-forward  /  backpropagation  neural  network  (NN),  were  employed  to  join 
together  the  molecular-scale,  particle-scale,  and  continuum  models  begun  in  Phase  I.  In 
addition,  a  large-scale  organic  gel  database  (over  39  different  entries)  was  created  to  train 
and  validate  the  NN.  To  facilitate  operator  handling,  the  network  mid  database  were 
wrapped  in  a  convenient  and  robust  Graphical  User  Interface  (GUI),  written  in  the 
platform  independent  Python  programming  language.  Preliminary  results  show  that  this 
software  (named  SOGeF-Synthesis  and  Optimization  of  Gel  Formulations)  offers  a 
reliable  and  comprehensive  design  tool  for  the  prediction  of  gel  formulations.  SOGeF 
uses  gel  constituent  properties  (i.e.  solvent  density,  boiling  point,  silica  content,  etc.)  to 
make  predictions  of  gelled  properties.  This  approach  is  advantageous  over  physics-based 
models  in  that  input  parameters  are  easily  identified  and  can  be  arbitrarily  chosen.  In 
addition,  as  the  size  of  the  training  database  grows,  SOGeF  predictions  become  more 
reliable.  This  comprehensive  approach  accomplishes  the  Army’s  objective  of  a  model  to 
describe  the  physical  properties  of  gelled  propellants. 
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Preface 


This  is  the  final  report  for  the  US  Army-RDECOM  ACQ  CTR  STTR  Phase  II  project 
entitled  “Computational  Design  Tool  for  the  Synthesis  and  Optimization  of  Gel 
Formulations  (SOGeF)”.  This  project  was  sponsored  by  the  US  Army-RDECOM  ACQ 
CTR  (Contract  Number  W91 1NF06C0186),  and  performed  by  researchers  at  CFD 
Research  Corporation  (CFDRC  Project  Number  8731).  The  principal  investigators  were 
Dr.  Stelu  Deaconu  mid  Dr.  Glen  Wilt  (CFDRC),  co-investigators  were  Dr.  Jeff  Morris 
(City  College  of  New  York)  and  Dr.  Boris  Khusid  (New  Jersey  Institute  of  Technology) 
and  the  project  manager  was  Matt  Thomas  (CFDRC).  The  Army-RDECOM  technical 
monitor  of  this  project  was  Dr.  Ralph  Anthenien  and  the  contract  representative  was 
Leroy  Hardy. 
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INTRODUCTION 


1  . 

1.1  Motivation  and  significance  of  the  current  research 

Gel  propulsion  systems  combine  the  best  characteristics  of  solid  and  liquid  propellants. 
The  gel  system  stores  like  solid  propellant,  but  flows  like  a  liquid  when  pressurized, 
enabling  throttle  and  restart  capability  similar  to  liquid  propellants.  In  addition,  gels  have 
a  lower  vapor  pressure  compared  to  their  liquid  counterparts.  Consequently,  the  storage 
containers  holding  them  can  be  much  lighter. 

Gelled  propellants  are  physical  gels  that  are  formed  and  stabilized  in  the  solvent  mainly 
through  physical  interactions  (secondary  forces  like  hydrogen  bonding,  van  der  Waals 
forces,  hydrophobic  forces,  etc.)  between  polymer  molecules  and/or  particles.  The 
common  feature  of  physical  gels  is  a  fragile  space- filling  connected  network.  Here  we 
focused  entirely  upon  gelation  caused  by  the  addition  of  solid  particles,  in  particular 
fumed  silica  (FS),  with  the  goal  of  being  able  to  predict  thermodynamic  and  rheological 
properties  of  this  network  based  on  knowledge  of  the  component  liquids  and  particles. 

The  physical  interactions,  which  cause  gelation,  will  impact  upon  the  preparation, 
loading,  storage,  flow,  injection,  atomization,  and  combustion  of  a  gelled  propellant  in  a 
propulsion  system.  The  properties  required  for  a  propellant  may  be  obtained  by 
manipulation  of  the  chemical  architecture  and  chemical  environment  of  its  components. 
In  Phase  I,  we  focused  primarily  on  the  gel  performance  in  storage  and  transportation, 
where  the  material  must  have  sufficient  rigidity  to  maintain  solids  in  suspension  (added 
metallic  powders,  for  example),  and  also  be  stable  against  network  collapse  and 
syneresis.  Our  Phase  I  emphasis  was  on  gelation  by  fumed  silica  for  a  range  of  liquid 
propellant  and  oxidizer  formulations.  The  propellant  and  oxidizer  systems  were  treated  as 
individual  gels,  i.e.  systems  where  mixing  occurs  only  as  part  of  the  combustion  process 
was  considered. 


Figure  1 :  View  of  Gelled  HMH/1RFNA  Combustion  Inside  a  Vortex  Motor  115  ms  After 
Ignition.  The  St  riot  ions  on  the  Upper  Lefi  Hand  Comer  of  the  Window  are  Deposits  of 
Silica.  The  Development  of  an  Innovative  Computational  Design  Tool  for  the  Synthesis 
and  Optimization  of  Gel  Formulations  (SOGeF)  in  Phase  II  has  made  Propellant 
Formulation  Development  More  Systematic. 
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An  enabling  technology  for  the  further  advancement  of  gelled  propulsion  is  the 
development  of  tools  to  render  the  synthesis  of  propellant  formulations  more  systematic. 
The  efforts  of  this  project  have  established  a  Quantitative  Structure  Property  Relationship 
(QSPR)  software  tool  (named  SOGeF-Synthesis  and  Optimization  of  Gel  Formulation) 
that  can  be  used  to  predict  gel  physical  properties  from  the  gel’s  constituent  properties. 

A  key  advantage  of  SOGeF  is  that,  apart  from  training,  model  predictions  are  explicit  and 
do  not  require  lengthy  physics-based  gelation  simulations.  In  addition,  SOGeF  input 
parameters  are  easily  identified,  as  apposed  to  the  challenging  specitication  physical 
constants  in  gelation  models.  Also,  as  the  size  of  the  training  database  grows,  predictions 
from  SOGeF  will  only  become  more  reliable,  wheras  the  accuracy  of  physics-based 
models  is  relatively  fixed. 

1.2  Proposed  solution:  OSPR-Neural  Network 

The  overall  flow  diagram  of  the  SOGeF  software  is  given  in  Figure  2.  SOGeF  operates 
by  puling  select  training  data  from  a  gel  property  database,  forming  intrinsic  relationships 
between  the  gel  descriptors  (constituents)  and  predictors  (gel  properties),  and  storing  the 
relationships  in  a  functional  form  (QSPRs)  for  future  predictions.  The  QSPR  selected  for 
this  study  was  that  of  a  feed-forward  /  backpropagation  neural  network. 

The  training  data  set  can  be  created  from  computed  and  /  or  experimental  data.  In  our 
Phase  II  implementation,  emphasis  was  placed  on  generating  experimental  data.  CFDRC 
teamed  up  with  the  New  Jersey  Institute  of  Technology  (NJIT)  and  City  College  of  New 
York  (CCNY)  to  measure  the  properties  of  over  39  different  organic-based  gels.  These 
were  used  to  train  SOGeF  in  two  distinct  categories  of  gels  (9  and  11%  silica  content  by 
weight). 


Figure  2:  SOGeF  Program  Flowchart 


2 


873 1/Final  Draft 


2  .  PHASE  I  TECHNICAL  OBJECTIVES  AND  ACCOMPLISHMENTS 

2.1  Phase  I  Technical  Objectives 

The  primary  technical  objective  of  the  Phase  I  program  was  the  development  of  a 
model(s)  to  describe  the  physical  properties  of  gelled  propellants.  To  obtain  these  goals 
requires  an  understanding  of  the  physico-chemical  structure  of  propellants  and  its  relation 
to  the  various  material  properties  that  takes  us  beyond  the  identification  of  mere 
empirical  trends. 

Specific  objectives  of  the  Phase  I  research  included: 

i.  Identification  of  critical  parameters  which  affect  the  structural  and 
material  property  characteristics  of  gelled  propellants 

ii.  Provide  the  basis  for  modeling  gelled  propellant  performance  based  on  the 
component  properties 

2.2  Summary  of  Phase  I  Accomplishments 

Experimental  studies  were  conducted  at  CCNY  on  the  rheological  properties  of  the  gels 
while  NJIT  conducted  syneresis  studies.  Given  the  difficulties,  both  infrastructural  and 
logistic,  associated  with  mixing  and  shipping  samples  of  hypergolic  gels  to  the  various 
research  institutions,  it  was  decided  to  work  with  gel  simulants  in  Phase  I.  The  simulant 
chosen  for  these  preliminary  studies  was  a  gel  composed  on  fumed  silica  and  water.  To 
ensure  consistency  and  statistical  reproducibility,  all  samples  were  prepared  by  CFDRC 
in  single  batches  and  subsequently  divided  and  sent  to  the  various  labs  for  testing.  All 
gels  contained  Cabot  Cab-O-Sil  M5  fumed  silica.  Four  different  types  of  gel  were 
prepared.  The  first  was  a  standard  8:1  water  silica  gel  (8  parts  of  water  to  1  part  of  fumed 
silica  by  weight).  CFDRC  also  mixed  8.25:1,  8.5:1,  and  9:1  silica  gels  so  that  the  effects 
of  silica  concentration  could  be  evaluated.  These  are  denoted  Gel  8,  Gel  8.25,  Gel  8.5, 
and  Gel  9,  respectively. 

2.3  Rheology  of  Fumed  Silica  Gels 

Rheological  studies  on  samples  of  aqueous  fumed  silica  gels  prepared  with  Cabot  CAB- 
O-SIL  M-5  silica  were  carried  out.  A  complete  battery  of  rheometric  tests  of  the  fumed 
silica-water  gels  was  conducted.  The  rheology  of  the  gels  showed  an  apparent  yield 
stress  of  several  hundred  Pascals  and  an  elastic  modulus  of  the  order  of  104  Pa;  the 
breakdown  strain  of  these  materials,  as  well  as  the  weaker  gels  prepared  onsite,  is  of  the 
order  of  0.06. 

Gel  samples  were  also  prepared  at  City  College  for  the  purpose  of  estimating  the  gel 
point  of  the  fumed  silica-water  system.  The  two  sets  of  samples  are  not  directly 
comparable,  however;  a  sample  with  a  water-to-silica  ratio  of  9:1  prepared  by  this 
procedure  formed  a  weaker  gel  than  the  9:1  sample  (Gel  9)  prepared  by  CFDRC. 
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Some  experiments  were  carried  out  in  a  controlled-stress  rheometer  (TA  AR  1000)  using 
40  mm  diameter  parallel  plates.  As  discussed  below,  the  results  show  that  the  samples 
prepared  by  CFDRC  are  strong  gels  with  elastic  moduli  of  about  10,000  Pa  and  yield 
stresses  of  about  450  Pa  (both  values  depending  on  concentration).  The  16.6:1  water-to- 
silica  gel  seems  to  be  very  close  to  the  point  of  incipient  gelation. 

Figure  4.1  shows  the  storage  modulus  (G’)  and  loss  modulus  (G”)  for  Gel  9  as  functions 
of  frequency  in  the  controlled  stress  rheometer.  LVR  data  was  taken  after  application  of  a 
slow  steady  shear  of  0.02  /  s  for  two  and  five  minutes  (total  strains  of  2.4  and  6)  to  Gel  9. 
The  data  for  both  moduli,  G’  and  G”,  essentially  superimpose  after  this  pre-shear  as  seen 
in  Figure  3. 


Figure  3:  G’  (filled  symbols)  and G”  (open  symbols)  for  Gel  9 following  the  same 
loading  procedure  with  an  additional  pre-shear  at  0. 02  /s  for  2  and  5  minutes.  The 
response  is  characteristic  of  a  strong  gel,  with  the  elastic  component  (G  ’)  independent  of 
frequency  and  more  than  an  order  of  magnitude  larger  than  the  viscous  component  (G”). 

The  response  is  characteristic  of  a  strong  gel,  with  the  elastic  component  (G’) 
independent  of  frequency  and  more  than  an  order  of  magnitude  larger  than  the  viscous 
component  (G”).  Strain  was  monitored  and  plotted  as  a  function  of  time.  Results  are 
shown  for  Gel  8.5  in  Figure  4. 
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Figure  4:  Step- stress  and  recovery  curves  at  different  imposed  shear  stresses,  Gel  8.5. 

Both  G’  and  G”  decrease  dramatically  at  higher  strains,  which  is  certainly  due  to 
breakdown  of  the  gel  structure  ( 

Figure  5).  Similar  behavior  is  observed  in  Gels  8.25  and  8.  The  first  deviation  in  G’  from 
the  low-strain  limit  occurs  at  a  strain  of  about  4%  for  Gel  8.5;  the  modulus  G’  is  11,000 
Pa,  corresponding  to  a  yield  stress  of  440  Pa.  The  high-  and  low-frequency  asymptotes  of 
G’  intersect  at  a  strain  of  6.6%,  which  would  give  a  yield  stress  of  720  Pa.  The  maximum 
in  G”  occurs  at  a  strain  of  17%,  corresponding  to  a  yield  stress  of  1500  Pa  based  on  the 
low-frequency  asymptote  of  G’  (G’  and  G”  are  nearly  equal  at  this  strain,  with  a  value  of 
340  Pa;  the  magnitude  of  the  complex  modulus  G*  =  V  (G’2  +  G”2)  is  470  Pa.).  The 
corresponding  values  of  the  yield  stress  for  Gel  9  are  deviation  from  linearity,  260  Pa; 
intersection  of  asymptotes,  1500  Pa;  and  the  maximum  in  G”,  1700  Pa. 


Figure  5:  Linear  moduli  versus  strain  at  a  frequency  of  1Hz  for  Gel  9  (triangles)  and  Gel 

8  (squares). 
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The  stress  relaxation  modulus  G(%  t)  following  step  strains  of  0.5,  1,  and  5%  on  the 
ARES  rheometer  is  shown  in  Figure  6  for  Gel  9.  The  data  overlap  for  y=  0.5%  and  1%, 
establishing  that  1%  is  within  the  linear  range  but  5%  is  not.  This  indicates  that  there  are 
short-time  processes  acting  on  a  separate  length  scale  that  affect  the  onset  of  nonlinearity; 
these  are  likely  to  involve  rapid  deformation  of  the  equilibrium  structure. 


t(s) 

Figure  6:  Relaxation  modulus  of  Gel  9  as  a  function  of  time  following  step  strains  of  0.5 

(circles),  1  (triangles)  and  5%  (squares). 

The  dynamic  moduli  of  Gel  IB,  a  14.7:1  water-to-silica  ratio  gel  prepared  at  City 
College,  are  shown  in  Figure  7  at  1  Hz  as  functions  of  strain.  The  moduli  are  several 
orders  of  magnitude  smaller  than  moduli  for  the  9:1  gel  prepared  by  CFDRC  ( 

Figure  5).  This  is  because  the  data  presented  are  for  lower  solids  content,  as  we  sought  to 
determine  the  gel  point.  Interestingly,  the  onset  of  nonlinearity  occurs  at  roughly 
comparable  strain  levels  in  the  weak  gels  and  strong  gels,  suggesting  that  strain  and  not 
stress  may  be  the  critical  variable  defining  failure. 
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Figure  7:  Storage  (circles)  and  loss  moduli  (triangles)  at  1  rad/s  for  Gel  IB  versus  strain. 


G’  and  G”  for  Gel  1C,  16.6:1  water -to-silica  ratio,  are  plotted  together  in  Figure  8.  The 
moduli  follow  nearly  parallel  power -law  behavior  over  the  two  decades  of  accessible 
frequency  range.  According  to  the  criterion  of  Winter  and  Chambon  (1986),  this 
concentration  is  close  to  the  incipient  gel  point. 


Figure  8:  Storage  and  Loss  Moduli  of  Gel  1C  As  Functions  of  Frequency. 


2.4  Svneresis  of  Fumed  Silica  Gels 

Experiments  were  conducted  on  water-based  silica  gels  provided  by  CFDRC.  A  DuPont 
Sorvall  RC  28S  SUPRAspeed  Centrifuge  was  used  to  measure  the  syneresis  of  water- 
based  silica  gels  GEL  8,  GEL  8.25  and  GEL  8.5  provided  by  CFDRC  over  the  range  of 
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accelerations  from  500  G  to  30,000  G.  Experiments  were  conducted  at  20°C.  As  an 
example,  photos  presented  in  Figure  9  show  a  GEL  8  sample  before  and  after  the  30-min 
exposure  to  5,500  G.  We  weighted  the  initial  sample  and  the  exudated  water.  Data  on 
syneresis  were  then  expressed  in  terms  of  the  mass  of  exudated  water  as  the  fraction  of 
the  initial  sample  mass. 


Figure  9:  Photos  of  a  GEL  8  sample  before  ( left j  and  after  (right)  the  90-m  in  exposure  to 
5,500  G  at20oC.  The  circle  indicates  the  location  of  water  exudated  from  the  gel. 

Plots  in  Figure  10  show  the  dependence  of  the  amount  of  exudated  water  as  the  fraction 
of  the  initial  sample  mass,  y,  on  the  magnitude  of  the  acceleration  following  the  30-min 
exposure.  The  extrapolation  of  logarithmic  fitted  curves  to  y  =  0  indicates  that  the 
syneresis  would  start  at  the  acceleration  of  about  2226  G  for  GEL  8,  590  G  for  GEL  8.25, 
and  526  G  for  GEL  8. 5. 

Plots  in  Figure  11  show  the  syneresis  of  GEL  8  as  a  function  of  the  exposure  time  for  the 
accelerations  5,500  G,  15,500  G,  and  23,500  G.  For  comparison,  we  also  measured  the 
time  dependence  of  syneresis  of  GEL  8.25  and  GEL  8.5  for  the  acceleration  of  5,500G  by 
varying  the  exposure  time  from  1  min  to  30  min.  The  mass  fraction  of  water  exudated 
from  GEL  8  during  the  90-min  exposure  to  1,500  G  was  found  to  be  about  0.28%.  Plots 
in  Figure  12  present  data  on  syneresis  of  GEL  8,  GEL  8.25  and  GEL  8.5  for  the 
acceleration  of  5,500G. 
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%  Syneresis  (mass)  vs  Gravity 


Figure  10:  Syneresis  of silica  gels  GEL  8,  GEL  8.25  and  GEL  8.5  for  200C following  the 
30-min  exposure  to  acceleration  (in  G).  The  amount  of  exudated  water  is  given  as  the 

fraction  of  the  initial  sample  mass. 


%  Syneiesis  (man)  vs  Time 


Time  (mini 

Figure  11:  Time  dependence  of  syneresis  of  GEL  8  following  the  exposure 

to  acceleration  for  20oC. 
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%  Synei  esis  (mass)  vs  Time 


Time  {mint 

Figure  12:  Time  dependence  of syneresis  of  GEL  S,  GEL  8.25  and  GEL  8.5  following  the 
exposure  to  acceleration  of 5,5000  for  20oC. 

To  visualize  the  network  formed  by  silica  particles  in  a  water-based  silica  gel?  we 
removed  water  using  the  supercritical  extraction  process  in  liquid  carbon  dioxide  [Taylor 
(1996)].  The  water  and  C02  are  fully  miscible  under  these  conditions,  so  that  the  water- 
C02  surface  tension  is  zero.  Thus  capillary  forces  that  would  collapse  the  pores  and 
rearrange  the  particles  during  the  water  extraction  are  avoided. 

Experiments  were  conducted  on  samples  of  gels  GEL  8,  GEL  8.25  and  GEL  8.5.  Table  1 
summarizes  the  data.  The  reproducibility  of  measurements  was  studied  by  repeating 
experiments  (denoted  as  Run  1  and  2  in  Table  1)  and  quantifying  several  SEM  images 
taken  from  the  same  sample. 


Gel 

Sample  No 

Void  shape 
factor 

Mean  void 
diameter,  imt 

Area  fraction 
of  voids 

GEL  8 

Run  1 

0.66 

86.29 

0.434 

GEL  8.25 

Run  1 

0.71 

98.72 

0.592 

0.64 

92.86 

0.625 

Run  2 

0.68 

95.32 

0.601 

GEL  8.5 

Run  1 

0.54 

171.22 

0.617 

0.62 

154.22 

0.585 

Run  2 

0.41 

150.10 

0.641 

GELS 

30-min  exposed  to 
11,500  G  at  20° C; 
syneresis 

33.12fv/v)% 

Run  1 

0.68 

74.25 

0.144 

0.69 

80.88 

0.176 

Run  2 

0.64 

63.46 

0.132 

0.68 

73.80 

0.159 

Table  1:  Characteristics  of  structures  formed  by  silica  particles 
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The  SEM  image  in  Figure  1 3  shows  the  network  structure  formed  by  silica  particles  in 
GEL  8.  The  surface  fraction  of  pores  and  the  mean  diameter  of  pores  are  found  to  be 
about  43%  and  86  ran,  respectively  (Table  1.1).  The  SEM  images  presented  in  Figures 
14  and  15  show  structures  formed  by  silica  particles  in  samples  GEL  8.25  and  GEL  8.5 
revealed  by  the  supercritical  extraction  of  water.  The  surface  fractions  and  the  mean 
diameter  of  pores  in  these  samples  are  found  to  be  about  59-62%  and  95-99  ran  for  GEL 
8.25  and  58-64%  and  150-171  nm  for  GEL  8.5  (Table  1.1).  This  explains  why  these 
samples  being  subjected  to  centrifuging  lose  water  faster  than  GEL  8.  (Figures.  10  and 
12) 

The  bulk  density  of  the  structure  formed  by  silica  particles  after  the  supercritical 
extraction  of  water  from  the  inline  filter  loaded  with  a  GEL  8  sample  was  found  to  be 
527.76  kg/m3.  Also,  we  attempted  to  characterize  the  rearrangement  of  silica  particles 
caused  by  centrifuging  a  gel.  The  mass  of  water  exudated  under  these  conditions 
constituted  about  16.2%  of  the  initial  mass  of  the  sample.  The  SEM  images  presented  in 
Figure  13  reveal  structures  formed  by  the  silica  particles  as-received  and  exposed  to 
centrifuging. 


The  SEM  image  of  the  structure 
formed  by  particles  in  GEL  8 
revealed  by  supercritical 
extraction  of  water  in  liquid 
carbon  dioxide. 


The  SEM  image  of  the  structure 
formed  by  particles  in  GEL  8.25 
revealed  by  supercritical 
extraction  of  water  in  liquid 
carbon  dioxide. 


The  SEM  image  of  the  structure 
formed  by  particles  in  GEL  8.5 
revealed  by  supercritical 
extraction  of  water  in  liquid 
carbon  dioxide. 


Figure  13 


Figure  14 


Figure  15 


11 


8731/FmalDraft 


Figure  16:  The  SEM  images  of  structures  formed  by  particles  in  GEL  8  revealed  by 
supercritical  extraction  of  water:  (left)  as-received;  (right)  following  the  30-min  exposure 
to  11,500  G  at  200C,  syneresis  33. 12  (v/v)  %. 


We  employed  nuclear  magnetic  resonance  (NMR)  spectroscopy  to  probe  the  interaction 
of  water  molecules  with  silica  particles  [Clifford,  1968;  Zimmerman,  1957],  Water 
molecules  are  dynamically  exchanged  between  an  environment  in  which  they  relax 
slowly  (free  water)  and  one  in  which  they  relax  rapidly  (bound  water).  Under  conditions 
of  rapid  exchange  between  free  water  and  bound  water  on  the  NMR  time  scale,  a  single 
relaxation  rate  is  observed,  which  represents  the  weighted  average  of  the  individual  rate 
for  each  environment.  We  used  17 O  nuclear  magnetic  relaxation  to  probe  the  mobility  of 
water  molecules  in  silica  gels.  For  comparison.  Oxygen- 17  NMR  measurements  were 
carried  out  for  an  empty  glass  tube,  a  tube  filled  with  D2O  up  to  2  inches  of  height,  a  tube 
loaded  with  D2O  and  8  wt/wt  %  silica  powder,  and  a  tube  loaded  with  D2O  mixed  with 
GEL  8.  The  chemical  shift  (relative  delta  in  resonance  frequency  reported  in  ppm)1  of  the 
170  resonance  peak  appears  to  be  about  3.5  ppm  for  D2O,  3.7  ppm  for  D2O  with  silica 
powder,  and  4.7  ppm  for  D2O  with  GEL  8,  as  shown  in  Figure  17  These  preliminary 
results  indicate  that  NMR  relaxation  can  be  employed  to  characterize  the  interactions  of 
solvent  molecules  with  particles  in  silica-based  gels. 


1  8  —  (v  -  Vrf.f)  xl(f  /  Vrf.f 
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Figure  17:  Oxygen- 17  NMR  resonance  peaks  for  (1)  D20,  (2)  D20  with  DegussaA200 

silica  powder,  and  (3)  D20  with  GEL  8. 

2.5  Molecular  Level  Modeling  and  Simulations 

As  experimental  data  were  available  for  adsorption  of  water  on  silica,  this  system  was 
used  as  a  test  to  validate  quantum  mechanical  and  simulation  methods  that  will  be  used  in 
Phase  II  for  propellant  (IRFNA,  MM  A  and  DMAZ)  +  fumed  silica  systems.  The 
procedure  to  be  used  for  those  systems  is  the  one  described  here  for  the  aqueous  filmed 
silica  systems.  The  force  field  for  aqueous  fumed  silica  gels  was  developed  and  validated 
using  GCMC  (Grand  Canonical  Monte  Carlo)  simulation  to  predict  water  adsorption  on  a 
silica  nanoparticle.  The  tasks  completed  in  Phase  I  included:  (1)  a  search  of  the  literature 
for  the  suitable  interaction  potentials  to  describe  Si02  and  water  interactions;  (2)  an 
evaluation  of  the  of  the  Si02-Si02  interaction  potential  by  Pereira  et  al.(1999)  and  van 
Beest  et  al.(1990)  that  have  been  published;  and  (3)  development  of  an  Si02  +  water 
interaction  potential  as  this  was  not  available  in  the  literature.  The  interaction  energy 
between  atoms  in  the  Si02  force  fields  is  of  the  form: 

r  r 

V  V 

where  Uff  is  the  potential  energy  between  atoms  i  and  j.  A,  b  and  C  are  parameters,  q  is 
the  charge.  We  use  this  potential  for  each  atom  in  Si02.  We  have  used  the  extended 
simple  point  charge  (SPCE)  model  [Berendsen,  1987]  for  water,  as  this  model  accurately 
represents  the  phase  envelope  (density  versus  temperature)  of  water  when  used  in 
simulation. 

We  developed  a  force  field  for  Si02  +  water  (cross  interaction)  using  computational 
quantum  mechanics.  The  interaction  between  the  solid  Si02  and  H20  was  calculated 
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using  the  model  clusters  for  Si02,  Si(OH)4,  Si20(0H)6,  Si303(0H)6,  and  SuO^OFQs  as 
shown  in  Figure  1 8. 


Figure  18:  Four  different  structures  of  silicon  clusters.  The  red  atoms  are  terminal 
oxygen  atoms  (only  bound  to  1  Si  atom),  the  blue  atoms  are  bridging  oxygen  atoms 
(bound  to  two  Si  atoms)  and  the  grey  atoms  are  silicon  atoms. 

The  electrostatic  potential  energy  surface  for  the  above  isolated  clusters  was  obtained 
from  quantum  mechanical  calculations.  The  charges  on  each  of  the  atoms  of  Si02  were 
determined  by  fitting  point  charges  at  each  atomic  center  to  this  electrostatic  potential 
energy  surface.  The  interaction  energies  of  each  of  these  Si02  clusters  and  a  water 
molecule  were  then  calculated  as  a  function  of  the  distance  between  the  centers  of  mass 
of  the  molecules  at  different  relative  orientations  of  the  Si02  cluster  and  the  water 
molecule.  The  interaction  energies  obtained  are  shown  in  Figure  19. 


Figure  19:  Calculated  energies  of  the  various  molecular  arrangements  ofSi02  cluster 
and  water  as  a  function  of  center-  of -mass  separation  distance 

The  calculated  quantum  mechanical  energies  were  then  used  to  parameterize  the  force 
field  for  the  interaction  of  water  and  various  Si02  clusters  used  in  simulation. 

Experimental  data  on  the  adsorption  of  water  on  silica  powders  (fumed  silica)  were 
measured  by  Meyer  and  Hackerman  (1 966).  We  used  this  data  to  test  the  force  field  that 
we  have  developed  by  comparison  with  results  from  Grand  Canonical  Monte  Carlo 
simulations  of  the  adsorption  of  water  on  a  4  nm  silica  nanoparticle  in  an  8  nm  simulation 
box.  Figure  20  contains  snapshots  of  such  simulations  at  low  and  high  relative  humidity. 
The  model  particle  we  have  used  in  simulation  was  obtained  by  computer  annealing  of  a 
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spherical  particle  in  a  molecular  dynamics  simulation  with  the  force  field  by  Van  Beest 
(1990). 


Figure  20:  Snapshot  of  of  fined  configurations  a  silica  particle  with  water  adsorbed  to  its 
surface  at  low  and  hign  relative  humidities,  respectively. 

Figure  21  is  a  comparison  of  the  adsorption  isotherm  from  simulations  using  the  force 
field  we  developed  together  with  the  available  experimental  data  at  25°C.  The 
simulations  were  run  a  two  different  relative  humidities  P/Po  (where  Po  is  the  vapor 
pressure  of  pure  water).  The  agreement  with  experiment  is  very  good. 


Figure  2 1:  Adsorption  isothermfor  the  adsorption  of  water  onto  the  surface  of  a  silica 
particle  at  300  K  compared  to  the  experimental  data  [Meyer,  1 966]. 

2.6  Particle  Scale  Models  and  Simulation 


We  have  established  a  capability  to  simulate  strongly  attractive  spherical  particles  in 
viscous  liquids.  Using  heuristic  but  reasonable  criteria  based  on  closeness  of  particle 
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pairs  to  be  termed  “connected,”  we  have  determined  conditions  resulting  in  fully 
connected  structures,  a  condition  needed  for  support  of  a  shearing  force  without  flow. 
We  employ  simulations  that  may  incorporate  hydrodynamic  interactions  (HI)  through  a 
Newtonian  liquid,  and  that  account  for  inter-particle  forces  expressible  as  a  gradient  of  a 
pair-particle  interaction  potential  (and  thus  dependent  on  knowledge  of  particle  positions 
only). 

Note  that  accounting  for  HI  implies  that  the  motion  of  one  particle  influences 
neighboring  particle  motions  (and  is  a  many-body  problem  dependent  on  instantaneous 
positions),  whereas  a  simple  drag  that  is  independent  of  interaction  hydrodynamics  is  a 
single-body  problem.  Based  on  the  solid  understanding  of  this  condition  for  attractive 
particles,  we  have  focused  on  development  of  the  Stokesian  Dynamics  algorithm, 
implemented  in  standard  and  accelerated  (Sierou  and  Brady  2001)  form,  for  strongly 
attractive  conditions.  We  are  not  aware  of  prior  use  of  the  method  for  particles  with 
strongly  attractive  spherically  symmetric  potentials. 

All  of  the  potentials  considered  have  a  continuous  force,  unlike  the  square-well  potential 
that  is  used  in  some  studies  (Foffi  et  al.  2005).  This  is  illustrated  for  two  forms  of 
potential  in  Figure  22. 


r 


Figure  22:  Examples  of  interparticle  force  forms  employed,  plotted  as  a  function  of  r, 
with  r =1  being  particle  contact.  The  squares  denote  a  potential  with  a  generalized 
Lennard-Jones  (of  36-18  form)  soft-core  potential  and  minimum  ofV/kT  =  -9.5  at  r  = 
1.07,  while  the  crosses  show  a  combined  steeply  repulsive  hard-core  potential  near  r=l 

with  a  linear  attraction  out  to  r=1.2. 

Successful  simulation  with  these  and  other  potentials  has  been  demonstrated.  In  order  to 
probe  the  microscale  basis  of  the  rheology  of  the  gelled  state  in  attractive  colloids,  we 
have  developed  forms  with  well-defined  separation  of  inflection  in  the  attractive  portion 
from  the  energy  minimum,  and  forms  such  as  the  linear  attraction  in  which  the  inflection 
can  be  made  arbitrarily  close  to  the  minimum. 
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We  find  that  strong  attraction,  with  V  >  5  kT,  results  in  structural  development  similar  to 
that  reported  in  Foffi  et  al  (2005)  for  Brownian  dynamics  simulations.  The  energy 
function  used  is  the  combined  form  discussed  above  and  illustrated  in  Figure  23  with 
minimum  of  9.5  kT.  In  simulations  in  three  dimensions,  we  have  gone  to  1000  particles. 


Figure  23:  Plot  of  particle  positions  after  100  Brownian  diffusion  times  (100a2 /Do) 
showing  clustering  in  a  15%  by  area  dispersion  with  strong  attraction. 

An  issue  that  was  proposed  for  study  was  the  role  of  shear  in  setting  structure,  and  we 
have  demonstrated  using  small  systems  (less  than  100  particles)  at  larger  fractions  that 
the  local  structures  generated  can  be  affected  by  even  a  relatively  weak  shear,  see  Figure 
24. 


Figure  24:  Full  symbols  show  structure  in  unit  cell  of  a  small  two-dimensional  attractive 
particle  simulation,  with  minimum  energy  of  -3  kT  atr/a  =  2. 05,  forPeclet  number  ofPe 
=  0.1  (Pe  =  ya2 /Do  where  y  is  the  shear  rate)  after  a  strain  of  10. 
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3  .  PHASE  II  TECHNICAL  OBJECTIVES  AND  ACCOMPLISHMENTS 

3.1  Objectives  of  the  Phase  II  study 

The  overall  objective  of  this  project  was  to  develop  Quantitative  Structure-Property 
Relationships  (QSPR)  for  the  prediction,  synthesis,  and  optimization  of  gelled  propellant 
formulations.  The  specific  objectives  of  the  Phase  II  effort  were: 

1 .  Prepare  a  number  of  gels  of  various  solvent  types  and  silica  contents 

2.  Create  an  extensive  gel  property  database,  including  the  following  measured 
properties  of  1 : 

a.  Density  (CFDRC) 

b.  Modulus  of  elasticity  (CFDRC) 

c.  50-caliber  test  (CFDRC) 

d.  Syneresis  (NJIT) 

e.  Yield  stress  and  loss  /  storage  modulus  (CCNY) 

3.  Develop  a  QSPR  model  (neural  network)  to  predict  future  gel  properties  from 
those  already  observed  in  the  database 

4.  Create  a  Graphical  User  Interface  (GUI)  for  the  QSPR  model  to  facilitate 
algorithm  usage  and  database  management 

5.  Continue  efforts,  from  Phase  I,  to  model  the  structural  formation  of  gels  through 
Brownian  /  Stokesian  dynamics 

6.  Study  the  dynamics  of  gelation  through  Dynamic  Light  Scattering  (DLS)  and  time 
dependent  rheology 

3.2  Summary  of  the  Phase  II  accomplishments 

The  objectives  of  the  Phase  II  study  were  successfully  accomplished.  These 
accomplishments  can  be  summarized  as  follows: 

1.  Gel  mixing 

Mixing  techniques  were  established  for  the  preparation  of  organic-based  gels. 
Over  39  different  gels  were  then  prepared  and  tested.  Consistency  error  estimates 
were  determined  for  various  preparation  uncertainties,  which  were  shown  to  be 
negligible  (below  measurement  uncertainty). 

2.  Gel  property  database 

Measurement  techniques  were  established  and  perfected  over  the  coarse  of  the 
Phase  II  project.  Gel  properties  were  measured  to  form  a  training  and  validation 
database  for  the  QSPR  model.  The  database  now  contains  23  different  organic 
entries  from  gels  created  with  9%  silica  and  16  entries  from  gels  created  with  1 1% 
silica. 

3.  OSPR-Neural  Network  Development 

A  feed-forward  /  back  propagation  neural  network  was  created  to  predict 
properties  of  new  gel  formulations  from  those  observed  in  the  gel  database.  The 
algorithm  was  demonstrated  with  excellent  results  over  the  relatively  limited 
training  dataset  (in  terms  of  typical  neural  network  training).  Suggestions  were 
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offered  and  demonstrated  to  deal  with  issues  arising  from  limited  amounts  of 
training  data,  which  show  substantial  possibilities  to  improve  network 
predictability. 

4.  Graphical  User  Interface  (GUI)  Development 

A  GUI  was  developed  for  the  QSPR  model  in  Python  programming  language, 
which  is  platform  independent.  The  GUI  facilitates  algorithm  usage  and  database 
management. 

5.  Brownian/ Stokesian  gelation  dynamics  model 

The  inter-particle  force  potential  was  established  and  demonstrated  for  aqueous 
colloidal  silica  dispersions.  A  procedure  for  simulating  larger  systems  faster  and 
without  much  loss  in  accuracy  was  selected  (11). 

6.  Dynamic  gelation  measurements 

Both  DLS  and  rotational  rheometry  were  used  to  study  the  formation  of  aqueous 
gels.  DLS  provided  insight  into  the  structural  evolution  of  gels,  while  rotational 
rheometry  explored  their  mechanical  evolution.  These  measurements  were  used 
to  determine  the  onset  of  gelation  and  compared  quite  well. 

Phase  II  efforts  have  successfully  produced  a  validated  QSPR  (neural  network)  model 
capable  of  making  accurate  gel  predictions  in  a  user  friendly  way.  Details  of  the  work 
performed  during  the  Phase  II  project  are  presented  in  detail  in  the  following  sections. 
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4. 


TECHNICAL  APPROACH 


4.1  Measurement  Procedures 

The  workload  of  preparing  and  testing  the  gels  was  divided  between  each  of  the 
participating  contractors.  The  following  sections  describe  the  experimental  procedures 
that  were  utilized.  Each  organization’s  responsibilities  are  outlined  below: 

□  CFDRC 

o  Gel  Mixing 
o  Density 

o  Modulus  of  Elasticity 
o  50-caliber  Test 

□  NJIT 

o  Syneresis 

o  Extrusion  and  Drop  Formation 
o  Evaporation 

□  CCNY 

o  Rotational  Rheometry 
o  Gelation  Dynamics 

4.1.1  Density7 

Gel  densities  were  measured  with  a  10.49-ml  Fisher  Scientific  pycnometer  and  a  Mettler- 
Toledo  Analytical  Balance  (see  Figure  25).  The  total  estimated  uncertainty  (instrument  + 
operator  variability)  in  measured  density  was  less  than  ±1%.  In  addition,  measurements 
were  made  by  CFDRC  in  Huntsville  AL,  where  the  samples  were  prepared. 
Consequently,  the  gels  did  not  incur  significant  degradation  due  to  long  storage  times. 


Figure  25:  (A)  Fisher  Scientific  Pycnometer,  and  (B)  Digital  Precision  Balance  Used  to 

Determine  Gel  Density 


4.1.2  50  Caliber  Test 

The  50-caliber  test  was  based  on  specifications  sent  by  Scott  Michaels  (AMRDEC).  A 
gel  fails  the  test  if  it  will  flow  through  a  14-inch  hole  with  a  head  height  of  7  inches; 
otherwise  the  gel  passes.  The  intent  of  the  test  is  to  simulate  a  ruptured  Hellfire  sized 
missile  lying  on  its  side. 
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Initially,  CFDRC  prepared  gels  with  9%  silica  content  by  weight.  Nearly  all  of  these  gels 
failed  the  50-caliber  test.  The  silica  content  was  then  raised  to  11%,  resulting  in  a 
substantial  pass-rate  improvement.  Table  2  summarizes  the  results  of  the  50-caliber  test. 

Table  2:  Summary  of  50-Caliber  Test  Results  Showing  a  Much  Improved  Pass  Rate  in  the 

11%  Data  Over  the  9%  Data 


Solvent 

9%  Silica 

iso-Pentane 

Fail 

Heptane 

Fail 

Octane 

Fail 

iso-Octane 

Fail 

CycloHexane 

Fail 

Hexene 

Fail 

Benzene 

Fail 

Diethyl  Amine 

Pass 

n-PropylAmine 

Fail 

Cyclohexene 

Fail 

11%  Silica 

n-Pentane 

Fail 

Pass 

Hexane 

Fail 

Pass 

Decane 

Fail 

Pass 

MEK 

Fail 

Fail 

Acetone 

Fail 

Fail 

MIBK 

Fail 

Fail 

Decanol 

Fail 

Fail 

Methanol 

Fail 

Fail 

NitroMethane 

Pass 

Pass 

NitroProp 

Fail 

Pass 

NitroEthane 

Fail 

Pass 

Triethyl  Amine 

Pass 

Pass 

1 -Propanol 

Fail 

Fail 

Dibutyl  Amine 

Fail 

Fail 

1-Pentanol 

Fail 

Fail 

n-ButylAmine 

Fail 

Fail 

4.1.3  Modulus  of  Elasticity 

The  modulus  of  elasticity  of  solids  and  semi-solid  materials  can  be  calculated  with  the 
following  formula 


E  =  c2-p 

Where  c  is  the  speed  of  sound  in  the  material  and  p  is  its  density.  For  our  gels  the  density 
was  available  from  measurement  (as  described  in  Section  4. 1. 1).  The  speed  of  sound  was 
found  by  measuring  the  time  required  for  an  ultrasonic  pulse  to  travel  a  known  distance 
through  the  gel.  The  apparatus  to  determine  the  speed  of  sound  is  shown  in  Figure  26. 
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Figure  26:  Setup  for  Gel  Sound  Speed  Measurements  (1)  DP  2000  Ultrasonic  Doppler 
Velocimeter  /  Pulse  Generator,  (2)  Tektronics  TDS-2 10  Digital  Oscilloscope,  (3)  Dial 
Caliper,  (4)  Gel  Samples,  (5)  Ultrasonic  Probe,  and  (6)  Plexiglass  Disk 


During  these  measurements  the  gel  was  placed  in  a  thin  polyethylene  bag  (small 
sandwich  bag)  and  set  between  the  tip  of  the  ultrasonic  probe  and  a  Plexiglass  disc 
(precision  machined  for  a  smooth  flat  surfaces).  An  ultrasonic  8-MHz  signal,  emitted  by 
a  5-mm  diameter  cylindrical  probe,  was  coupled  (acoustically)  to  the  gel  bag  by  a  thin 
film  of  GE  Ultrasonic  Gel  (also  used  to  couple  the  bag  to  the  Plexiglass  disk).  The  speed 
of  sound  (c)  was  then  calculated  with: 


_  2  d_ 

At 

Where  d  is  the  distance  between  the  tip  of  the  ultrasonic  probe  to  the  Plexiglass  disk  and 
At  is  the  time  between  the  outgoing  ultrasonic  pulse  and  the  reflection  from  the  first 
Plexiglass  surface.  The  sound  speed  uncertainty  was  on  average  ±4%  of  the  recorded 
value,  which  includes  error  introduced  by  the  wall  thickness  of  the  bag  and  operator 
errors  made  during  reading  the  oscilloscope  trace.  In  addition,  the  sound  speed 
measurement  of  each  gel  sample  included  10-12  repeated  measurements  at  distances 
ranging  from  0.2  to  0.7  in. 

4.1.4  Syneresis 

Samples  of  organic-based  fumed  silica  gels  were  sent  by  CFDRC  to  the  NJIT  for 
syneresis  testing  (the  complete  set  of  syneresis  data  is  located  in  Appendix  B).  A  DuPont 
Sorvall  RC  28@  SUPRAspeed  centrifuge  was  used  to  induce  syneresis  in  the  gels  over 
the  range  of  accelerations  from  500  to  30,000  G  at  20  °C.  Figure  27  shows  an  example 
of  hexane  silica  gel  both  before  and  after  a  30-min  exposure  to  10,000  G.  As  can  be  seen 
in  the  image  on  the  right,  solvent  has  been  exudated  from  the  sample  through  the  collapse 
of  the  gel  structure  (i.e.  reduced  gel  volume). 
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Figure  27:  Photos  of  Hexane  Silica  Gel  Before  (left)  and  After  (right)  a  30-min  Exposure 
to  10,  OOOG  at  20°C.  The  Circle  Indicates  the  Location  of  Organic  Solvent  Exudate  d 

From  the  Sample. 

Syneresis  was  expressed  as  the  ratio  of  exudated  solvent  mass  to  initial  sample  mass. 
Measurement  repeatability  was  established  by  duplicating  experiments  under  typical 
operating  conditions.  Table  3  summarizes  the  consistency  in  duplicated  measurements. 
It  is  noted  that  the  greatest  variation  between  experiments  occurs  at  the  smaller 
acceleration  (500  G),  due  to  a  small  value  of  exudated  mass.  It  is  also  noted,  that  these 
values  are  not  related  to  repeatability  over  long  sample  storage  times,  which  will  be 
addressed  at  the  end  of  the  Section. 


Table  3;  Reproducibility  of  Syneresis  Measurements 


Solvent 

Operating  Conditions 

Syneresis,  mass% 

Variation, 

% 

Acceleration 

(G) 

Temp 

(°C) 

Time 

(min) 

Run  1 

Run  2 

Heptane 

20,000 

20 

10 

36 

37 

3.0 

Decane 

500 

20 

90 

0.7 

0.8 

9.4 

Decane 

30,000 

20 

5 

37 

39 

-3.7 

In  Phase  I,  the  syneresis  of  water-based  silica  gels  was  found  to  be  a  strong  function  of 
the  exposure  time  to  acceleration.  Conversely,  the  syneresis  of  organic-based  gels  did 
not  demonstrate  a  strong  dependence  on  exposure  time  at  fixed  accelerations  (see  Table  4 
and  Figure  28).  Once  subjected  to  centrifuging,  an  organic  gel  looses  a  certain  amount  of 
solvent  very  fast,  but  then  retains  the  remaining  amount  of  solvent. 

Notice  that  some  values  in  Table  4  show  a  re- absorption  of  the  solvent  (within  the 
inherent  measurement  error).  This  may  indicate  a  time  dependent  swelling  of  the 
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collapsed  gel  structure  after  the  removal  of  acceleration  (similar  to  a  sponge). 
Consequently,  residual  swelling  should  be  more  pronounced  under  smaller  accelerations, 
where  a  greater  percentage  of  elasticity  is  recovered.  It  may  be  necessary  to  study  this 
phenomenon  further  in  order  to  reduce  uncertainty. 


Table  4:  Syneresis  of  organic  solvents  gelled  with  9  w/w%  fumed  silica  following  the 
ex£osure  to  50CJG,  2, 000G  and  30,000G  accelerationsjor  5  min  and  90min  at  20°C. 


Gels  (9% 
w/w  silica) 

Syneresis,  mass  % 

500  G 

2,000  G 

30,000  G 

5  min 

90  min 

5  min 

90  min 

5  min 

90  min 

Iso-Pentane 

0.9 

0.7 

0.6 

1.1 

35.4 

38.5 

n-Pentane 

0.3 

0.4 

0.4 

0.6 

30.0 

29.1 

Hexane 

0.2 

0.1 

0.5 

3.0 

37.8 

40.0 

Heptane 

0.1 

0.1 

3.3 

5.2 

35.7 

40.0 

Octane 

0.1 

0.1 

6.8 

8.4 

38.2 

42.6 

Iso-Octane 

1.7 

1.5 

12.4 

11.5 

41.6 

46.7 

Decane 

0.7 

0.7 

9.7 

9.9 

36.9 

42.8 

%  syneresis  (mass)  vs  Time  (5,000G) 
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Figure  28:  Time  dependence  of  syneresis  in  organic  solvents  gelled  with  9  w/w%  fumed 
silica  subjected  to  5,000G  acceleration  at  20°C. 


Similar  to  the  water-based  gels  of  Phase  I  (see  Figure  30),  the  syneresis  of  organic-based 
gels  showed  strong  dependence  on  acceleration  rate.  Figure  29  shows  the  syneresis  of 
the  alkanes  following  5  and  10-minute  exposures  to  various  acceleration  rates.  The 
syneresis  tends  to  taper  off  at  high  acceleration  rates,  indicating  that  the  gel  is 
approaching  a  limit  of  compression  that  cannot  be  recovered. 
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%  syneresis  (mass)  vs  Gravity  (5-min) 
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Figure  29:  Syneresis  of  organic  solvents  (iso-pentane,  n-pentane,  hexane,  heptane, 
octane,  iso-octane,  and  decane)  gelled  with  9  w/w%  fumed  silica,  following  the  5  (top) 
and  10  (bottom)  minute  exposure  to  acceleration  (in  G)  at  20°C. 
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%  Syneresis  (mass)  vs  Gravity 


Figure  30:  Syneresis  of silica  gels  SILICA  50 ,  SILICA  8.25  (8.25 parts  of  water  to  1  part 
of  fumed  silica  by  weight)  and  SILICA  8.50  ( 8.50 parts  of  water  to  1  part  of  fumed  silica 
by  weight)  following  the  30-min  exposure  to  acceleration  (in  G)  at  20°C.  The  amount  of 
exudated  water  is  given  as  the  fraction  of  the  initial  sample  mass.  [7] 

To  reduce  the  50-caliber  failure  rate,  CFDRC  increased  the  silica  loading  of  the  gels  from 
9  to  11%.  The  effect  of  this  action  on  syneresis  can  be  seen  in  Figure  31.  It  appears  that 
syneresis  decreases  with  increasing  silica  content.  Also,  MIBK  demonstrated  the  highest 
reduction  in  syneresis. 
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Figure  31:  Syneresis  of  organic-based  silica  gels  with  9  and  11%  (w/w)  silica  following 
the  exposure  to  acceleration  20,000  (top)  and  10,000  (bottom)  Gfor  30min  at  2(fC. 


The  data  collected  in  Figure  31  was  extracted  from  relatively  fresh  batches  of  gels  on 
different  test  dates.  Consequently,  storage  time  had  little  effect  on  the  results  shown. 
The  adverse  effect  of  long  sample  storage  times  on  syneresis  repeatability  is  shown  in 
Table  5.  A  reduction  in  syneresis  with  time  is  attributed  to  evaporation.  Under 
evaporation,  the  collapsed  gel  volume  will  be  similar,  but  the  exposed  solvent  volume 
will  not. 
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Table  5:  Comparison  of  Syneresis  after  Long  Storage  Time  (Decane) 


Gels  (9%  silica  w/w) 

%  Syneresis  (mass)  of  2,000G  and  90min 

Decane 

August  2007 

March  2008 

9.86 

6.35 

4.1.5  Extrusion  and  Drop  Formation 

The  jetting  characteristics  of  CS  gels  were  studied  with  the  setup  shown  in  Figure  32. 
Syringes  1  and  2  (both  10-ml  with  16 -mm  barrel  IDs)  are  coupled  together  by  Tygon 
tubing,  filled  with  oil  to  monitor  the  pressure  acting  on  the  piston  rod  of  syringe  2.  The 
piston  rods  of  syringes  2  and  3  are  directly  connected  together.  Organic  gel  was  loaded 
into  syringe  3  (glass  5 -ml  with  11. 6 -mm  barrel  ID).  It  was  necessary  to  use  a  glass 
syringe  to  hold  the  gel,  since  rubber  piston  seals  would  swell  and  eventually  cause  the 
piston  to  stall.  Contrary,  a  glass  syringe  does  not  react  with  the  organic-based  gel  and 
provides  a  relatively  low  friction  between  the  piston  rod  and  syringe  barrel. 


Syringe  2 


Syringe  3 


Syringe  1 


Pressure  Gauge 
(0-100  psi) 


Clamping  system 
aligned  vertically 


Harvard  Apparatus 
PHD2000  Syringe 
Pump  for  generating 
pressure. 


Figure  32:  Gel  Jetting  Setup:  (Left)  Experimental  setup  for  testing  pressure  driven 
syneresis  (Right)  Syringe  3  loaded  with  gel.  (Volume:  5ml;  barrel  ID:  11.6mm,  nozzle 

ID:  2.0  mm,  nozzle  length:  17.0  mm) 


The  extrusion  tests  included  two  steps.  First,  a  blank  test  (empty  syringe  3)  was 
conducted  at  the  desired  piston  speed.  This  was  done  to  measure  the  driving  pressure  due 
to  friction  in  the  setup  (oil  and  piston  seals).  Next,  syringe  3  was  loaded  with  a  gel  and 
the  test  was  conducted  at  the  same  piston  speed  as  the  blank  case.  The  difference  in 
pressure  between  the  loaded  and  blank  case  was  recorded  as  the  driving  pressure  of  the 
gel. 


The  setup  was  first  used  to  study  the  flow  of  a  hexane-based  silica  gel.  Following  the 
application  of  ~1 5  psi  a  solvent  drop  appeared  from  the  nozzle  (Figure  33).  Solvent  and 
CS  particles  slowly  filled  the  pendant  drop  causing  it  to  eventually  detach  from  the 
nozzle.  An  increase  of  pressure  (to  ~26  psi)  caused  a  transition  from  dripping  to  a 
column  extrusion. 
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Figure  33:  Pressure  Driven  Flow  of  a  Hexane-based  Gel  Showing  Drop  Formation  at  15 

psi 

Figure  34  shows  the  dripping  (left)  and  jetting  (right)  action  of  a  CS  gel  (4.7%)  flowing 
through  two  different  orifices  at  400  ml/min.  With  a  2 50 -micron  diameter  orifice  (left), 
the  gel  drips  out.  With  a  2 -mm  orifice  at  the  same  flow  rate  (right),  the  gel  extrudes  as  a 
column.  Notice  the  light  regions  in  the  figure,  indicating  the  presence  of  silica. 


Figure  34:  Extrusion  and  Drop  Formation  of  a  Colloidal  Silica  Gel  (4.7%)  at  the  Same 
Flow  Rate  Through  a  250-micron  (Left)  and  2-mm  (Right)  Orifice 
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4.1.6  Gel  Structure  During  Solvent  Evaporation 

Preliminary  studies  on  the  evolution  of  the  gel  structure  during  solvent  evaporation  were 
conducted  on  the  NJIT  fluorescence  microscope  system  shown  in  Figure  35.  The 
NIKON  microscope  is  equipped  with  two  fluorescence  CCD  cameras  (a  high-resolution 
CoolSNAP  SQ2  and  a  high-speed  Cascade  128+),  a  dual-view  aperture  imaging  system, 
and  a  fluorescence  illumination  system  EXFO  X-Cite  120.  The  system  allows  the 
simultaneous  recording  of  two  types  of  fluorescent  markers. 


Two  CCD  Cameras: 
Photometries 
CoolSNAP  SQ2 
(shown  in  picture) 
and  Cascade  128+, 
which  can  collect 
more  than  500  full 
frames  of  1 6-bit 
data  per  second. 


Air  Table 


Computer/ 

Data  Requisition 


Figure  35:  NIKON  fluorescence  microscope  system  used  to 


Dual-View  Aperture 
Imaging  System 

Microscope  (Nikon 
SMZ1 500)  with  5  Filters 

Fluorescence 
Illumination  System 
(EXFO  X-Cite  120) 

Stage  with 
Position  Monitor 


Figure  36  shows  the  structure  transformations  of  a  hexane-based  silica  gel  during  the 
evaporation  of  hexane.  As  the  hexane  vaporized,  the  gel  structure  began  to  shrink  and 
break  up.  Shrinkage  of  the  silica  network  occurs  when  the  imbibitional  pressure  created 
by  the  solvent  is  removed.  The  gel  structure  transformations  can  be  quantified  by 
embedding  fluorescent  markers  in  the  gel. 
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4.1.7  Rheometrv  (Yield  Stress.  Storage,  and  Loss  Modulus) 

All  rotational  rheometry  was  performed  with  a  Rheometrics  ARES  strain-controlled 
rheometer  at  25  °C.  The  rheometer  works  by  measuring  the  response  (shear  stress:  x)  of 
the  gel  to  an  oscillatory  strain  (s).  In  purely  elastic  materials,  x  occurs  in  phase  with  s.  In 
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purely  viscous  materials,  x  will  lead  s  by  90°.  The  relationship  of  shear  stress  to  strain  of 
an  intermediate  material  (elastic  +  viscous)  is  given  by: 

s  =  s0  sin(<y  •  t ) 
r  =  t0  sin(<y  -  t  +  S) 

From  this,  the  shear  storage  (G’)  and  loss  (G”)  modulus  are  defined  as: 

Storage  :  G'=  —  cos(e>) 

Loss  :G"=  —sin(S) 

A  sweep  of  frequency  was  performed  at  each  strain  amplitude  to  determine  G’  and  G”. 
The  limit  of  the  gel’s  linear  response  (yield  point)  was  acquired  by  increasing  strain 
amplitude  until  the  shear  stress  began  to  disappear.  At  this  point,  yield  stress  (  r  =  G'-s  ) 

was  determined.  If  the  strain  did  not  exceed  the  established  yield  point,  the  test  could  be 
performed  again  with  the  same  plateau  in  x;  otherwise  the  plateau  would  always  be  lower 
(confirming  the  established  assessment  of  yield  point). 

The  values  of  G’,  s,  and  x  (at  the  yield  point)  were  reported  for  each  of  the  organic  gels 
tested  and  formed  the  basis  of  subsequent  property  prediction.  Figure  37  shows  the 
response  of  a  Tri-ethylamine  gel  (11%  Cabosil)  to  increasing  strain  amplitude.  This  gel 
had  a  yield  stress  of -460  Pa,  as  seen  in  the  peak  of  the  s  •  G'  curve. 
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Strain  (%) 


Figure  3  7:  Storage  Modulus  G  Loss  Modulus  G  ’  and  Stress  vs.  Strain.  Yield  Stress 

~460Pa.  (T ri-ethylamine  11%  Cabosil) 


The  primary  determinant  of  gel  yield  point  properties  was  found  to  be  silica-loading  level 
(see  Figure  45).  Yield  stress  increased  with  silica  fraction  ( c ),  within  the  alkanes, 
approximately  as: 


zy=K(c-coy 

Where  c0  denotes  the  minimum  silica  weight  fraction  capable  of  producing  observable 
gel  properties.  For  the  alkanes:  a  &  1  and  K  was  on  the  order  of  25,000-30,000  Pa.  Not 
all  gels  had  the  same  measurable  characteristics. 

Strong  gels  of  varying  strength  were  formed  by  the  alkanes,  alkenes,  and  cyclic  alkanes. 
Small  nitro-genated  alkanes  (nitro-methane,  nitro-ethane,  and  nitro-propane)  all  formed 
strong  gels,  with  nitro-methane  being  notably  stronger  than  the  others  (+50%).  Weaker 
gels  were  formed  from  acetone  (order  of  magnitude  weaker  than  alkanes).  Additionally, 
all  of  the  alcohol-based  gels  were  found  to  be  weak  and  failed  to  remain  gelled  unless 
stored  with  no  ullage  space,  to  prevent  evaporation. 

The  general  conclusion  from  the  rheometry  studies  is  that  the  presence  of  an  -OH 
functionality  in  an  organic  solvent  results  in  poor  gelation  (see  Figure  38).  The  fact  that 
ketone  with  =0  (doubly  bonded  oxygen)  can  undergo  the  keto-enol  tautomerization  (8)  to 
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form  an  -OH,  under  the  influence  of  a  trace  acid  or  base,  suggests  that  the  surface  of  the 
fumed  silica  may  interact  with  ketones  to  promote  this  reaction.  The  evidence  to  this 
point  seems  to  support  literature  claims  (9)  that  the  stabilizing  force  is  due  to  solvation 
effects  associated  with  the  tendency  of  materials  with  H-bonding  capabilities  (other  than 
water)  to  coat  the  silica  surface. 


o  o  o  o 


SOL  (stable) 


Silica  particles 


Particles  interact 


(B)  Weakly  H-bonding  liquid 


by  H-bondlng 

oo 

13- 


gel 


Figure  38:  Schematic  representation  of  two  possible  scenarios  that  can  occur  in  the  case 
of  silica  particles  dispersed  in  a  liquid:  (A)  strongly  H-bonding  liquid  (a  solvation  layer 
is  formed  on  the  particles  by  H-bonding,  resulting  in  a  stable  sol);  (B)  weakly  H-bonding 
liquid  (the  particles  interact  directly  by  H-bonding  to  form  a  gel)  [9J. 
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4.1.8  Gelation  Dynamics 

The  gelation  process  was  studied  on  aqueous  colloidal  silica  (CS)  dispersions  comprised 
of  Ludox  SM  30  (Grace-Davison  /  Sigma- Aldrich).  The  spherical  CS  particles,  in  these 
dispersions,  are  in  the  15-18  nm  diameter  range,  based  on  light  scattering  analysis  in 
diluted  samples.  Dispersions  were  received  as  a  stable  mixture  of  30%  silica  (by  weight) 
in  water.  This  corresponds  to  a  baseline  silica  volume  fraction  (<|))  of  0. 16;  the  fraction 
was  either  reduced  or  increased  through  the  addition  of  water  or  controlled  evaporation 
respectively. 

As  delivered,  the  CS  dispersions  have  a  pH  of  9  or  slightly  greater.  At  this  pH,  the  silica 
surface  is  negatively  charged,  stabilizing  the  dispersion  against  flocculation  (providing  a 
protective  maximum  in  the  interaction  potential  due  to  Coulombic  repulsion).  The 
attractive  forces,  which  bind  the  particles  of  a  gel  structure  together,  are  thought  to  be  due 
to  van  der  Waals  interactions  (which  occur  beyond  the  Coulombic  barrier);  but  may  also 
involve  hydrogen  bonding  of  the  silica  surfaces.  To  destabilize  the  dispersion  in  a 
controlled  and  tunable  fashion,  varying  levels  of  aqueous  NaCl  solution  was  added  to  the 
Ludox.  Brinker  and  Scherer  (3)  indicate  that  the  added  ionic  strength  of  the  solvent  upon 
salt  addition  causes  the  Coulombic  repulsion  from  ionized  surfaces  of  the  dispersed  silica 
particles  to  be  screened,  leading  to  capture  in  the  attractive  potential  (van  der  Waals). 
The  observed  behaviors  (structural-light  scattering  and  mechanical-rheology)  following 
the  mixing  time  and  quantity  of  added  NaCl  solution  is  described  further  in  the  following 
sections. 

4. 1.8.1  Dynamic  Light  Scattering  (DLS) 

Dynamic  Light  Scattering  (DLS)  was  used  to  access  the  structure  (average  flocculated 
object  size)  of  the  emerging  gel.  DLS  measures  the  de-correlation  with  time  of  the 
diffuse  intensity  in  light  passing  through  the  mixture.  De-correlation  in  intensity  occurs 
near  the  mean  free  time  of  the  mixture,  which  can  be  linked  to  the  average  particle  radius 
through  Brownian  diffusion  relationships. 

The  DLS  studies  were  performed  using  scattering  at  90  degrees  from  the  incident  Ar  laser 
beam  (wavelength  of  488  nm)  in  either  fixed  or  mobile  mode;  the  latter  involving 
translation  at  speeds  on  the  order  of  microns/second  to  study  the  gelled  materials.  A 
standard  photo-multiplier,  counter,  and  correlator  system  was  used. 

Mixes  of  aqueous  CS  were  prepared  with  varying  levels  of  NaCl  solution  (used  to 
destabilize  the  dispersion)  and  monitored  by  DLS  for  average  hydrodynamic  radius 
(<r>).  As  can  be  seen  in  Figures  39  and  40,  <r>  has  an  initial  power  law  growth  period 
followed  by  a  rapid  divergent  growth  (indicating  the  onset  of  gelation).  In  addition,  the 
rapid  growth  rate  occurs  earlier  with  higher  NaCl  concentration,  due  to  a  breakdown  in 
the  Coulombic  repulsion  barrier.  This  effect  becomes  saturated  at  large  NaCl 
concentrations,  where  the  growth  is  limited  solely  by  diffusion.  It  is  also  noted  that  the 
mixtures  become  more  poly- disperse  as  the  flocculation  process  occurs  (up  to  the  point 
where  the  gel  network  begins  to  span  the  entire  mixture)  and  there  should  be  significant 
spread  around  the  mean  radius. 
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Figure  39:  Hydrodynamic  Radii  of  Clusters,  Measured  with  DLS,  During  the  Gelation 
Process  Following  Addition  of  Salt  (%  weight)  at  Various  CS  Percentages 


Figure  40:  Hydrodynamic  Radii  of  Clusters,  Measured  with  DLS,  During  the  Gelation 
Process  Following  Addition  of  Salt  (molarity)  at  Various  CS  Percentages 
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The  time  at  which  divergence  of  <r>  is  observed  is  one  measure  of  the  gelation  time  of 
the  dispersion.  Table  6  shows  DLS-based  measurements  of  gelation  time  for  some  of  the 
NaCl  concentrations  tested. 


Table  6:  DLS-Based  Measurements  of  Gelation  Time  at  Various  NaCl  Concentrations 
_ (From  Results  Shown  in  Figure  39) _ 


%  Weight  NaCl 

2.06 

1.97 

1.89 

1.78 

Gelation  Time  (min)  by  DLS 

40 

52 

92 

106 

A  second  more  fundamental  indication  of  gelation  time  is  obtained  through  dynamic 
measurements  of  elastic  modulus.  As  will  be  shown,  these  measurements  agree  well  with 
one  another,  although  the  rheological  approach  indicates  a  slightly  earlier  gelation  time. 
As  noted  in  Richter  et  al.  (4),  DLS-based  gelation  times  do  not  always  accord  with  the 
mechanical  definition  of  the  gel  determined  by  linear  viscoelastic  rheology.  This  is  likely 
related  to  the  fact  that  a  single  scattering  angle  experiment  probes  at  a  particular 
characteristic  length  (due  to  heterogeneity  of  the  gel)  while  rheometric  tests  intrinsically 
probe  all  scales.  Methods  (5)  have  been  applied  to  CS  gels  to  develop  improved 
understanding  of  the  heterogeneity  in  the  gel  structure,  which  could  help  establish  a  better 
approach  to  DLS-based  measurements  of  gelation  time. 

4. 1.8.2  Dynamic  Rheology 

The  same  materials  studied  by  DLS  have  been  monitored  continuously  using  oscillatory 
rheometry,  by  measuring  the  storage  (G’)  mid  loss  modulus  (G”).  In  Phase  I,  a  method 
using  an  overlaying  layer  of  low  viscosity  oil  (immiscible  with  the  aqueous  silica 
solutions)  was  shown  to  be  essential  for  extended  rheometric  measurements;  this 
approach  was  used  in  the  rheometry  described  here. 

The  tests  were  performed  at  a  fixed  strain  of  100%  per  cycle,  as  apposed  to  the  standard 
fixed  stress.  Figures  41  and  42  show  typical  responses  of  G’  and  G’  ’  after  the  addition  of 
salt.  The  values  of  both  G’  and  G”  are  seen  to  increase  strongly  at  a  time  starting  just 
under  80  minutes  after  salt  addition.  According  to  the  Winter  and  Chambon  (6)  criterion 
of  parallel  dependence  of  G’  and  G”  on  frequency,  the  gel  time  associated  with  each 
figure  appears  to  be  between  80  and  90  minutes.  The  corresponding  DLS  measured  gel 
time  of  Figure  39  is  noticeably  higher  (106  minutes). 


37 


8731/Final  Draft 


mm  M 


\m  l 


10  - 


j  - 


o  a  - 


o.oi 


o.  i 


T 

v 

* 

0 

4 

a 


Ci\  89  minutes 
ti".  89  minutes 
G\  90  minutes 
GM,  90  mmutes 
G\  105  minutes 
G\  105  minutes 


►  *0 


fi4A4»*tt****** 

*#•*  Saia*8o&e&S®£$ 

0oOOOOO°o0° 


•  • 


•  * 


A  a 

o 


7  V 

„  v£»T 


*#* 


r^ui 


10 


]00 


(o,  rads 

Figure  41:  The  storage  and  loss  moduli  obtained  at  a  strain  of  100%  from  the  dispersion 
of  1.78-wt%  salt  and  10.3  lwt%  Si02.  The  same  symbol,  open  and  closed,  respectively, 
is  used  for  G’  and  G”  at  the  noted  times  after  adding  the  destabilizing  salt  solution. 


iff* 


mini  - 


▼  G'  at  SO  minutes 
V  CT  at  SO  minutes 
*  G1  at  90  minutes 
o  G'  at  90  minutes 
■  ti1  at  145  minutes 
□  G1  at  145  minutes 


Cfl 

& 

b 


a 


UK)  - 


in 


n  □  □ 


* 

G 


□  □ 


□  □  □ 


□  □ 


□  □ 


^S8*fl8**n 


0000 


o  G  o  o 


OOOOO 


0.1  - 


;V^ 


VV, 


■v?¥ 


II. Ill 


m 


11111 


co  (rad/s) 

Figure  42:  Elastic  (G’)  and  loss  ( G ”)  moduli  a  a  Junction  of  the  oscillation  frequency  for 
a  CS  dispersion  of  5.2  %  silica  by  volume,  and  0.347  M  in  NaCl  at  three  times  following 

addition  of  salt. 
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It  is  noted  that  there  was  difficulty  in  establishing  the  exact  time  when  the  gel  criterion 
was  satisfied,  since  the  duration  of  the  frequency  scan  was  comparable  to  the  rate  at 
which  G’  and  G”  showed  significant  change.  In  addition,  the  LVR  rheological 
properties  were  at  the  limit  of  the  instruments  torque  measurement  accuracy  for  the  un¬ 
gelled  aqueous  CS  dispersion.  A  pragmatic  method  of  assessing  the  gelation  time,  which 
yields  an  essentially  equivalent  gel  time,  is  to  associate  gelation  with  the  onset  of  a 
growth  in  G’  to  reliable  measurement  levels. 

Following  gelation,  two  characteristics  of  the  gel  were  most  notable  (see  Figure  43).  The 
first  is  that  both  G’  and  G”  increase  markedly  with  time  following  the  onset  of  gelation. 
The  second  point  is  that  the  failure  strain  (s*)  decreases  considerably  (from  0.3  to  0.02 
%).  Loss  in  elasticity  of  the  gel  indicates  that  the  microstructure  is  increasingly  rigidified 
as  time  passes,  due  to  greater  connectivity,  which  allows  less  flexural  motion  without 
bond  rupture. 


— ■—  G’  190l 
□  G"190 

-L-G'  145( 
A  G"  14& 
frequency_1r 
vol_frac_4.70/ 
[NaCILO.347 


Figure  43:  Strain  Sweeps  at  190  and  1450  Minutes  After  Salt  Addition  ( 0.347  M)  for  a 

Ludox  Dispersion  at  <j)  =  0.047 
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4.2  Organic  Gel  Mixing  /  Database  Development 


The  main  goal  of  this  project  was  to  develop  a  QSPR  model  that  predicts  bulk  gel 
properties  from  the  various  components  and  mixing  techniques  used  in  creating  the  gel. 
Initial  efforts  were  aimed  at  linear  regression  relationships  originating  from  molecular 
level  descriptors,  such  as  molecular  electrostatic  potential  (MEP).  It  was  latter  decided 
that  using  indirect  macroscopic  descriptors,  such  as:  liquid  density,  viscosity,  and  silica 
loading,  would  work  equally  as  well  in  a  neural  network-based  architecture.  In  order  to 
train  the  network,  a  large  database  of  gel  properties  had  to  be  constructed. 

4.2.1  Initial  Production 

Organic  gels  were  prepared  by  CFDRC  in  rented  lab  space  at  Oakwood  College  and  latter 
at  Shelby  Research  Center  (Huntsville).  The  gels  were  mixed  (using  Cabosil  fumed 
Silica)  in  a  1 -gallon  pneumatically  actuated  mixer;  additional  smaller  batches  were  mixed 
with  a  500-ml  hand  actuated  mixer  (see  Figure  44). 


Figure  44:  (A)  Gallon  Size,  Pneumatic  Mixer,  and  (B)  500ml  Hand  Operated  Mixer 

(Images  Not  to  Scale) 

Initial  samples  of  gelled  heptane  were  prepared  to  form  a  baseline  for  subsequent  mixing 
techniques  and  Cabosil  loading  percentages.  With  the  product  literature  from  Cabot  and 
experience  in  making  water  gels,  an  initial  4%  Cabosil  loading  (by  weight  of  solvent) 
was  used.  This  resulted  in  a  gel  that  was  too  thin.  Cabosil  was  then  added  to  the  mixture 
(up  to  8.2%  total)  until  the  gel  was  too  thick  with  incomplete  mixing.  A  new  batch  of  gel 
was  prepared  with  7.6%  Cabosil,  which  was  again  too  thin.  Finally,  a  sequence  of  small 
batches  was  prepared  until  the  desired  consistency  was  attained  at  9%  Cabosil  loading. 

Following  the  heptane  loading  experiment,  CFDRC  concluded  that  after  a  gel  is  mixed,  it 
is  not  possible  to  disperse  any  more  Cabosil.  Consequently,  the  Cabosil  needs  to  be 
completely  wetted  with  the  liquid  solvent  quickly  (comparable  to  the  gelation  time)  to 
ensure  total  mixing  of  the  gel  batch. 
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As  with  heptane;  hexane,  pentane  mid  iso-pentane  were  all  gelled  with  9%  Cabosil. 
Gelling  of  alcohols,  however,  involved  more  than  simple  percentage  Cabosil  loading.  As 
mentioned  in  Section  3.1.7,  hydrogen  bonding  from  the  alcohol  prevents  the  development 
of  network  formations  between  Cabosil  aggregates.  This  can  be  overcome  through  the 
use  of  surfactants,  which  inhibit  bonding  between  the  solvent  and  Cabosil. 

The  surfactant  recommended  by  Cabot  (Duomeen  TDO)  is  a  waxy  solid  and  had  to  be 
completely  dissolved  into  the  solvent  before  mixing  in  the  Cabosil.  After  preparing  a 
number  of  small  test  batches,  it  was  found  that  an  8%  Cabosil  and  0.96%  surfactant 
loading  produced  the  desired  gel  consistency.  These  percentages  were  used  to  create  gels 
from  n-propanol,  iso-propanol,  methanol,  and  n-butanol. 

It  is  noted  that  the  total  required  loading  (surfactant  +  Cabosil)  of  alcohol-based  gels  was 
similar  to  the  alkanes  but  required  less  Cabosil.  Future  experimentation  may  find 
surfactants  capable  of  lowering  the  required  Cabosil  loading  even  further  for  all  types  of 
gels. 

4.2.2  Gel  Test  Matrix 

With  baseline  mixing  techniques  and  percentages  in  place,  an  initial  test  matrix  was 
established  to  determine  the  relative  influence  of  silica  loading  on  gels  of  the  following 
alkanes:  hexane,  octane,  and  decane.  The  gels  were  loaded  with  7,  9,  and  11%  fumed 
silica  mid  compared  by  yield  stress  using  data  in  the  Linear  Viscoelastic  Range  (LVR) 
(see  Figure  45).  It  was  noted  by  CCNY  that  the  gels  were  weak  and  failed  to  remain 
gelled  unless  stored  with  no  ullage  space. 


41 


8731/Final  Draft 


VJ 

*  ^ 
Vi 

TJ  S 

v  a 

Is  tfc 

S 

vS 

rjTi 

W 


1 

1  1  1 

Y  hexane 

□  octane 

1200 

x  decane 

n  - 

1000 

X 

000 

- 

600 

T 

□ 

400 

- 

?nn 

- 

n 

_ 1 _ 

,  <  , _ , _ 

0.02 


0.04 


0.06 


0.0S 


0.1 


0.12 


c  (wt  %  FS) 

Figure  45:  Results  of  Yield  Stress  far  Initial  Test  Matrix  of  Alkanes  (CCNY  data) 
thawing  a  Strong  Dependence  on  Silica  % 

As  can  be  seen  in  Figure  45,  silica  loading  is  more  influential  than  solvent  type  (likewise 
MEP  and  other  solvent  descriptors)  in  determining  yield  stress  within  the  alkanes  family. 
As  a  result,  it  was  determined  that  testing  only  select  members  of  each  organic  group 
would  significantly  reduce  the  size  of  the  test  matrix  without  loosing  generality  in  the 
capabilities  of  the  QSPR  model.  Consequently,  the  problem  of  determining  descriptors 
characteristic  of  a  specific  organic  group  would  have  to  be  resolved. 

In  order  to  draw  out  the  impact  of  solvent  type^  subsequent  gel  samples  were  mixed  in 
distinct  silica-loading  groups  across  all  of  the  selected  solvents.  The  first  group  that  was 
tested  (9%  silica)  resulted  in  gels  that  were  too  thin  (liquid  in  some  cases),  most  of  which 
failed  the  50 -caliber  test.  An  additional  group  (11%  silica)  was  added,  which  had  a  much 
higher  50-caliber  pass  rate  (except  for  the  alcohols,  ketones,  and  the  mono  and  di¬ 
sub  stituted  amines).  The  complete  gel  database  is  given  in  Tables  7  and  8.  Gels  from  the 
9%  data  that  were  too  thin  for  storage  modulus  measurements  were  assigned  a  default 
value  of  10  Pa  for  preliminary  QSPR  modeling  purposes. 
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Table  7:  The  (9%)  GelsDatabase  (06/08/2008) 


Measured  Gel  Properties  (9%  Silica  by  Weight)  1 

Storage  Modulus 

2000G/90min 

10000G/30min 

Density 

Bulk  Modulus 

.50  Cal  Test 

Solvent/Liquid 

G'  (Pa) 
CCNY 

(%)  Syneresis 

NJIT 

(%)Syneresis 

NJIT 

(g/ccm) 

CFDRC 

E  (MPa) 
CFDRC 

Pass/Fail 

CFDRC 

1 

N-Pentane 

36092 

0.6 

6.15 

0.71 

671 

Fail 

2 

Iso-Pentane 

21422 

1.1 

11.11 

0.69 

654 

Fail 

3 

Hexane 

20747 

3 

17.14 

0.76 

928 

Fail 

4 

Heptane 

17265 

5.2 

17.39 

0.76 

957 

Fail 

5 

Octane 

24501 

8.4 

21.69 

0.79 

1059 

Fail 

6 

Iso-Octane 

17270 

11.5 

25.75 

0.77 

883 

Fail 

7 

Decane 

19163 

9.9 

23.52 

0.83 

1189 

Fail 

8 

Cyclo-Hexane 

23400 

10.2 

30.58 

0.87 

1312 

Fail 

9 

Hexene 

14666 

8.57 

31.98 

0.77 

817 

Fail 

10 

Cyclo-Hexene 

32798 

2.21 

22.76 

0.91 

1288 

Fail 

11 

Benzene 

36014 

7.21 

27.48 

0.98 

1458 

Fail 

12 

Methanol  (w/Surf) 

729 

31.6 

46.72 

0.91 

1142 

Fail 

13 

1 -Propanol  (w/Surf) 

39 

39.1 

47.44 

0.91 

1367 

Fail 

14 

1 -Butanol  (w/Surf) 

10  (?)  Liq. 

38.34 

47.72 

0.93 

1370 

Fail 

15 

Hexanol  (w/Surf) 

10  (?)  Liq. 

41.11 

50.29 

0.94 

1528 

Fail 

16 

Octanol  (w/Surf) 

10  (?)  Liq. 

32.65 

39.55 

0.95 

1630 

Fail 

17 

Decanol 

4806 

12.96 

24.44 

0.93 

1662 

Fail 

18 

Nitromethane 

47170 

3.06 

23.91 

1.27 

2029 

Pass 

19 

Nitroethane 

29893 

7.68 

27.63 

1.16 

1723 

Fail 

20 

Nitropropane 

33383 

9.49 

28.98 

1.1 

1621 

Fail 

21 

Tri ethyl  Amine 

14300 

25.34 

43.5 

0.82 

1025 

Pass 

22 

Diethyl  Amine 

19000 

25.54 

43.15 

0.8 

1008 

Pass 

23 

n-PropylAmine 

10  (?)  Liq. 

41.04 

52.48 

0.82 

1248 

Fail 

24 

Acetone 

2357 

34.2 

43.46 

0.909 

1229 

Fail 

25 

Methyl-Ethyl  Ketone 

1600 

40.53 

50.5 

0.9 

1305 

Fail 

26 

MIBK 

655 

33.96 

63.18 

0.91 

1308 

Fail 

27 

1 -Propanol  No  Surf. 

1000 

39.1 

52.4 

0.91 

1324 

Fail 

28 

Octanol  No  Surf. 

N/A 

N/A 

N/A 

0.95 

1656 

Fail 

29 

1- Propanol  w.  Surf 

N/A 

N/A 

N/A 

0.92 

1346 

Fail 

30 

Butanol  No  Surf. 

N/A 

N/A 

N/A 

0.93 

1411 

Fail 

31 

Hexanol  No.  Surf. 

N/A 

N/A 

N/A 

0.94 

1540 

Fail 

Surfactant  -  Duomeen  0.9%  by  weight 


Table  8:  The  0 1%)  Gels  Database  (06/08/2008) 


Measured  Gel  Properties  (11%  Silica  by  Weight)  | 

Storage  Modulus 

2000G/90min 

10000G/30min 

Density 

Bulk  Modulus 

.50  Cal  Test 

Solvent/Liquid 

G1  (Pa) 
CCNY 

(%)  Syneresis 

NJIT 

(%)  Syneresis  (g/ccm) 
NJIT  CFDRC 

E  (MPa) 
CFDRC 

Pass/Fail 

CFDRC 

1 

N-Pentane 

58000 

0.12 

19.86 

0.72 

673 

Passed 

2 

Hexane 

36510 

7.17 

29.64 

0.77 

868 

Passed 

3 

Decane 

24000 

5.8 

22.96 

0.83 

1167 

Passed 

4 

MEK 

1689 

33.08 

43.17 

0.93 

1357 

Failed 

5 

Acetone 

3170 

27.22 

36.67 

0.92 

1258 

Failed 

6 

1-Pentanol 

310 

39.92 

49.56 

0.92 

1426 

Failed 

7 

MIBK 

894 

29.93 

41.26 

0.92 

1305 

Failed 

8 

Decanol 

3511 

4.07 

7.23 

0.95 

1731 

Failed 

9 

Methanol 

2089 

32.69 

42.12 

0.92 

1101 

Failed 

10 

Nitromethane 

74179 

0.03 

14.58 

1.29 

1964 

Passed 

11 

Nitropropane 

34887 

2.57 

20.26 

1.19 

1777 

Passed 

12 

Nitroethane 

62844 

0.03 

17.24 

1.19 

1751 

Passed 

13 

TriEthylAmine 

19813 

16.49 

37.2 

0.85 

1009 

Passed 

14 

n-ButylAmine 

1000 

34.25 

43.14 

0.84 

1289 

Failed 

15 

iso-PropylAlcohol 

5700 

30.88 

42.88 

0.9 

1176 

Failed 

16 

Di  butyl  Amine 

1000 

35.14 

45.63 

0.88 

1338 

Failed 
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4.2.3  Extraneous  Factors 

A  series  of  tests  were  conducted  to  determine  the  influence  of  extraneous  factors  on  gel 
properties.  These  include:  operator,  mixer  size,  surfactant  usage,  and  storage  time.  Of 
these  factors,  storage  time  (due  to  sample  transportation  and  lab-time  availability)  was 
found  to  be  the  most  influential  source  of  error. 

In  March  2008  the  NJIT  conducted  a  re-test  on  a  decane  sample  prepared  in  August  2007. 
The  results  showed  a  change  in  the  syneresis  (2,000G  /  90min)  from  9.9  to  6.35  %.  This 
change  was  attributed  to  the  loss  of  solvent  from  the  gel  due  to  evaporation. 

Repeat  samples,  not  adversely  influenced  by  storage  time,  were  prepared  for  the  solvents 
and  conditions  shown  in  Table  9.  The  average  density  (p)  and  modulus  of  elasticity  (E) 
of  the  hand  mixer  batches  were  found  to  be  within  ±  1%  of  the  respective  gallon  mixer 
batches.  This  inaccuracy  is  within  uncertainty  estimates  of  the  measurements  (±  1%-p,  ± 
5%-E).  Unfortunately,  the  volume  of  the  hand  mixer  batches  was  insufficient  to  be  sent 
to  CCNY  and  NJIT  for  testing.  However,  since  the  overall  appearance,  density,  and  E 
values  were  nearly  identical,  it  is  concluded  that  mixer  geometry,  mixer  size,  and 
variations  in  operator  induced  shear  rates  during  mixing  do  not  affect  gel  quality  and 
properties  significantly. 

As  mentioned  earlier,  surfactants  were  used  to  encourage  gelation  in  alcohol-based  gels. 
To  determine  the  impact  of  surfactant  usage,  mixes  of  short  (iso-propanol)  and  long 
(octanol)  chain  alcohols  were  prepared  with  and  without  surfactant,  using  both  the  gallon 
and  hand-mixer.  The  measured  values  of  density  and  E  were  identical  (within 
uncertainty)  for  each  case  of  solvent. 


Table  9:  List  of  Repeat  Mixes  Used  to  Determine  Extraneous  Factors  in  Gel  Preparation 


Solvent 

Gallon  Mixer  Batch 
(1) 

Hand  Mixer  Batch 
(2) 

Operator 
(1),  (2) 

n-Pentane 

G\  (%  S) 

p,  E 

WHD,  WHD 

Iso-Pentane 

G’,  (%  S) 

p,  E 

WHD,  SSD 

Hexane 

G\  (%  S),  p,  E 

p,  E 

WHD,  SSD 

Heptane 

G\  (%  S),  p,  E 

p,  E 

WHD,  WHD 

Octane 

G’,  (%  S) 

p,  E 

WHD,  WHD 

Iso-Octane 

G’,  (%  S) 

p,  E 

WHD,  SSD 

Decane 

G\  (%  S) 

p,  E 

WHD,  SSD 

Octanol  w/o  Surf. 

G’,  (%  S),  p,  E 

p,  E 

WHD,  WHD 

1-Propanol  w.  Surf. 

G\  (%  S),  p,  E 

p,  E 

WHD,  SSD 

(%  S)  -  (%)  Syneresis,  Surf  -  Surfactant  (here  0.9%  by  weight)  Duomeen  TDO  (N-tallow 
trimethylene  diamine  dioleate),  WHD  -  Waite  Dykes,  SSD  -  Stelu  Deaconu 
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4.3  OSPR  Model:  Neural  Network 


4.3.1  Theory 

The  established  method  for  estimating  properties  of  new  gel  formulations  is  based  on 
Quantitative  Structure-Property  Relationships  (QSPR).  Classic,  algebraic  QSPR  has 
been  successfully  used  in  the  past  for  pure  liquids.  Its  application  to  gels,  however,  was 
ambiguous  due  to  the  addition  of  extra  factors  (gellant,  surfactants,  etc.). 

The  chosen  method  was  to  implement  a  QSPR-like  property  prediction  algorithm  as  a 
feed-forward  /  backpropagation  neural  network  (NN).  A  literature  review  shows  that 
NN-based  QSPR  has  been  implemented,  with  good  results,  for  the  prediction  of  physical 
properties  of  organic  compounds  (1).  The  following  are  advantages  of  NN 
implementation: 

1)  Relatively  straightforward  computer  implementation 

2)  A  quantifiable,  unambiguous  breakdown  of  descriptor  influence  on  the 
structure-property  relationships 

3)  Use  of  solvent  macroscopic  physical  properties  (which  are  measured  rather 
than  computed) 

4)  The  possibility  of  comprehensive  learning  over  the  entire  liquid  /  solvent 
database  of  interest  (at  discrete  or  variable  gellant  mass  fractions) 

5)  Easy  to  train  /  re-train  with  various  (user)  data  sets. 

Figure  46  shows  the  overall  feed-forward  /  backpropagation  NN  architecture. 


Properties  From  Solvent  Descriptors 
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The  architecture  begins  with  values  of  select  properties  (descriptor  set)  passing  through  a 
scaler,  which  in  our  case  brings  values  into  the  (0,1)  range.  This  is  done  so  that  the 
network-training  algorithm  is  not  too  far  out  of  balance  upon  startup.  Next,  the  values 
are  weighted  (according  to  W)  mid  sent  to  the  various  nodes  of  the  hidden  layer.  The 
hidden  layer  nodes  accumulate  weighted  inputs  and  activate  outputs  ( z°lit )  through  an 

activation  function.  In  our  case,  the  Sigmoid  function  (shown  below)  is  used  for 
activation,  which  has  the  stabilizing  advantage  of  bounding  the  output  to  the  interval 
(0,1). 


Outputs  of  the  hidden  neurons  are  then  weighted  (according  to  V)  and  sent  to  the  nodes 
of  the  output  layer.  The  output  layer  nodes  accumulate  weighted  inputs  and  activate 
outputs  through  another  activation  function  (in  our  case  Sigmoid).  The  outputs  are  then 
re-scaled  (based  on  the  scale  utilized  during  training)  to  provide  the  predictor  set.  Output 
scaling  is  made  compatible  with  the  Sigmoid  function  range  (0,1).  In  addition,  the  scale 
is  set  in  excess  of  the  training  data  to  allow  for  extrapolated  predictions.  The  overall 
architecture  of  the  NN  can  be  described  by  the  following  equation: 


yf  =  f(Z  vu  f(Z  wjk  xk  ))>  ie  (t  °) 

7=1  k= 1 

Where  H,  I,  and  O  are  equal  to  the  number  of  hidden  neurons,  input  descriptors,  and 
output  descriptors  respectively. 

4.3.1. 1  Training 

The  main  problem  in  establishing  a  network  is  to  select  weight  arrays  (W  and  V),  which 
provide  mi  optimal  predictor  set  for  a  given  descriptor  set.  This  is  accomplished  by 
providing  positive  feedback  (training)  in  the  form  of  an  error  function.  The  error 
function  was  chosen  as  the  square  error  (SE)  between  the  predictor  ( ypred  )  and  training 

(ytrue )  set: 


SE  =  o.5^(yrd  ->T)2 


For  the  network  we  have  chosen,  the  sensitivity  of  SE  to  the  individual  weights  is  given 
by: 
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These  sensitivities  define  an  optimal  direction  for  weight  advancement  ( Aw„  and  Avy  ) 
given  by: 


A  dSE 

Aw  =  — 

8wIJ 

ass 

Ay,  =  -u - 

'7 

where,  ju  is  the  learning  rate.  In  our  implementation,  a  local  minimum  error  is  found 
along  the  optimal  direction  before  the  direction  is  updated  on  the  next  training  cycle. 
Training  proceeds  until  a  minimum  is  reached  in  the  SE  function.  However,  this 
minimum  is  not  guaranteed  to  be  optimal  over  the  entire  SE  function. 

4.3.1.2  Prediction  Risk:  Nonlinear  Cross  Validation  (NCV) 

As  mentioned,  the  error  function  is  used  to  establish  how  well  the  network  fits  its  training 
data.  A  prediction  risk  establishes  how  well  the  network  will  likely  predict  future  data. 
We  have  chosen  to  use  nonlinear  cross  validation  (NCV)  for  prediction  risk  estimation. 
In  NCV,  the  training  data  is  divided  into  v  disjoint  subsets  (D,  )  of  size  N, .  The 
network  is  trained  with  each  subset  removed  from  the  training  data,  which  forms  the 
following  error  (  CV  )  on  the  exempt  data: 


CK=±Z(y^-y“y 

A,  D , 

The  errors  of  all  exempt  subsets  are  then  averaged  to  yield  a  prediction  risk. 

4.3.1.3  Sequential  Network  Construction  (SNC) 

Another  important  issue  facing  network  selection  is  in  determining  its  size.  Excessively 
large  networks  tend  to  over-fit  the  training  data  at  the  expense  of  increased  prediction 
risk.  We  use  a  procedure  called  Sequential  Network  Construction  (SNC)  to  establish  the 
optimal  number  of  hidden  neurons  with  respect  to  NCV. 

In  SNC,  the  network  architecture  is  initialized  to  some  arbitrarily  small  number  of  hidden 
neurons.  After  initial  training  is  performed,  NCV  is  used  to  establish  a  prediction  risk. 
Nodes  are  then  added  to  the  hidden  layer  repeatedly,  with  training  originating  from  the 
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previously  determined  weights.  After  an  arbitrarily  large  number  of  hidden  neurons  has 
been  reached,  the  network  with  the  minimum  prediction  risk  is  chosen.  SNC  has  the 
following  main  advantages  (2): 

1)  As  larger  networks  are  constructed,  the  features  of  the  previous  networks  are 
preserved. 

2)  It  is  less  computationally  expensive  to  cultivate  a  large  network  from  a  small 
one  than  to  prune  nodes  of  a  larger  network  trained  from  scratch. 

The  coding  for  the  SNC  algorithm  used  in  this  study  is  listed  in  Appendix  C. 

4.3.2  Application 

4.3.2.1  Descriptors  /  Predictors 

For  the  NN  to  provide  worthy  predictions  beyond  curve  fitting,  there  must  be  intrinsic 
relationships  between  the  descriptors  and  output  variables.  As  demonstrated  earlier,  gel 
properties  follow  closely  with  silica  loading  and  less  with  the  solvent  properties  of  a 
given  organic  group.  Consequently,  appropriate  descriptors  should  be  characteristic  of 
the  organic  group  a  particular  gel  belongs  to,  with  some  variation  inside  each  group  to 
address  local  dissimilarities.  As  seen  in  the  distribution  of  MEP  based  parameters  for 
alkanes  and  alcohols  (Table  10),  some  MEP  variables  (i.e.  and  <7_2 )  would  work  well 

in  distinguishing  organic  groups.  However,  the  added  effort  of  computing  MEP 
descriptors  was  abandoned  for  the  more  computationally  economic  but  indirect 
descriptors  of  Table  11  (see  Appendix  A  for  9  and  11%  silica  loading  descriptors). 
These  were  established  by  manually  removing  low  impact  descriptors  from  a  larger 
arbitrary  property  matrix  in  order  to  optimize  NN  output  (more  than  50  descriptor 
combinations  were  tested).  NN  output  was  based  on  the  following  values: 


Storage  Modulus  G' 

(Pa) 

Yield  Stress 

(Pa) 

Syneresis  2000G/90min 

(%) 

Syneresis  10000G/30min 

(%)  ^ 

Gel  Density 

(gm/cm3) 

Modulus  of  Elasticity 

(MPa) 

Table  10:  MEP  Derived  Parameters  of  Alkanes  and  Alcohols  (Originally  Intended  as 
_ _ Descriptors.  in  the  QSPR  Model ')  _ _ 


Molecule 

■vs+ 

Vs“ 

<7_2 

_  2 

&1ot 

V 

Pentane 

8.1 

-5.0 

32.4 

10.9 

43.2 

8.1 

Heptane 

9.9 

-8.0 

28.8 

17.5 

46.3 

10.9 

Isopentane 

4.9 

-4.5 

13.6 

9.7 

23.4 

5.7 

Isopropyl 

alcohol 

6.1 

-11.6 

29.5 

114.0 

143.5 

23.4 

Methanol 

11.3 

-18.7 

28.8 

112.2 

141.1 

22.9 
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Butanol 

8.2 

-9.4 

39.1 

102.0 

141.1 

28.3 

Propanol 

5.6 

-11.3 

24.5 

106.0 

130.5 

19.9 

Table  1 1 :  Most  Effective  Descriptor  Set  (Solvent  Properties)  To  Date  (1/12/09)  Manually 


Selected_  to  Optimize  NN_  Output 


1.  Density 

(gm/ccm) 

2.  Viscosity 

(cPoise) 

3.  Surface  Tension 

(mN/m) 

4.  Dielectric  Constant 

5.  Vapor  Pressure 

(kPa  @  25°C) 

6.  Accentric  Factor 

7.  Dipole  Moment  (Debye) 

8.  Standard  Enthalpy  of  Formation  (HF) 

(kJ/mole) 

9.  Melting  Point 

(°C) 

10.  Boiling  Point 

_ m _ 

An  important  part  of  selecting  appropriate  descriptor  variables  was  to  determine  how  well 
the  resulting  NN  detected  patterns  in  the  data  (mainly  group  distinction).  Without 
significant  pattern  recognition,  the  NN  will  simply  curve  fit  the  training  data  and  exhibit 
limited  predictive  abilities.  Figure  47  shows  Percentage  Descriptor  Influence  (PDI) 
patterns  for  a  few  of  the  solvents  tested  (Alkanes  and  Alcohols  at  9%  silica  loading).  PDI 
was  calculated  by  individually  perturbing  each  input  variable,  of  a  given  solvent,  by  10% 
and  noting  the  percentage  change  in  the  output  array.  A  clear  pattern  can  be  seen  in  the 
PDFs  of  the  alcohols.  The  PDI  pattern  is  also  an  indication  of  a  group’s  sensitivity  to 
input  (predictor)  variable  error,  which  was  considered  negligible  during  this  study. 


49 


8731/Final  Draft 


(A)  N-Pentane 


(B)  Hexane 


(E)  n-Propanol 


(F)  Decanol 


Figure  47:  Percentage  Descriptor  Influence  for  Alkanes  (A,  B,  and  C)  and  Alcohols  (D, 
E,  and  F).  NN  Training  with  Set  of  Ten  Descriptors  and  Six  Gel  Parameters 
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The  PDI  patterns  of  Figure  47  can  be  an  indirect  indication  of  a  group’s  input  proximity. 
Solvents,  which  are  close,  will  display  similar  PDI  sensitivities.  Figure  48  shows  the 
average  Euclidian  distance  between  input  pairs  of  the  alkanes  and  alcohols.  As  can  be 
seen,  there  appears  only  a  slight  natural  grouping  of  the  descriptors  as  a  whole.  That  is, 
the  average  distance  between  any  two  input  descriptor  arrays  is  roughly  equivalent, 
regardless  of  the  two  groups  the  inputs  belong  to. 

Proximity  patterns  within  individual  input  variables  were  more  favorable  to  learning.  Of 
all  the  selected  inputs,  accentric  factor  was  found  the  most  capable  of  distinguishing 
groups  (see  Figure  49).  It  is  also  noted  that  alcohols  showed  the  greatest  PDI  with 
respect  to  accentric  factor. 


Figure  48:  Average  Euclidian  Distance  Between  Input  Pairs  of  the  Alkanes  and  Alcohol 
Groups  Showing  Only  a  Slight  Distinction  Between  the  Groups 
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Figure  49:  Average  Distance  (Accentric  Factor)  Between  Input  Pairs  of  the  Alkanes  and 
Alcohol  Groups  Showing  a  Significant  Distinction  Between  the  Groups 

4.3.2.2  Training  and  Validation 

Training  was  performed  under  constant  silica  loading  using  the  descriptor  (solvent)  and 
predictor  (gel)  sets  shown  in  Appendix  A  and  Tables  7  and  8  respectively.  A  NN  is 
typically  trained  with  a  large  set  of  well-distributed  data.  This  particular  study  was 
complex  in  that  the  distribution  and  number  of  training  data  are  both  limited  by  the  nature 
of  solvents  available  for  testing.  In  addition,  random  errors  within  the  training  data  can 
significantly  warp  the  results  of  the  network.  This  presents  difficulties  in  training  the  NN 
to  predict  rather  than  training  it  to  merely  curve  fit  training  sets.  In  spite  of  this,  the 
network  will  in  general  make  predictions  no  better  than  its  ability  to  fit  the  training  data. 
Figures  50  and  51  show  examples  of  the  NN’s  ability  to  fit  its  training  values. 
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Figure  50:  Comparison  ofNNFit  to  its  Training  Values  (Full  1 1%  Data  Base  with  23 
Hidden  Neurons)  Showing  Good  Agreement  Over  Most  of  the  Training  Data  (Red-Large 
Spread  in  the  Output  Data,  Yellow-Arbitrari ly  Chosen  Output  Values  Used  to  Fill  Void 

Data,  Blue-Normal  Data) 


Figure  51:  Comparison  ofNN  Fit  to  its  Training  Values  (Full  9%  Data  Base  with  23 
Hidden  Neurons)  Showing  Good  Agreement  Over  Most  of  the  Training  Data  (Red-Large 
Spread  in  the  Output  Data,  Blue-Normal  Data) 
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%  Error  in  Fit 


As  can  be  seen  in  Figures  50  and  51,  the  NN  is  capable  of  fitting  most  of  the  training  data 
very  well.  However,  some  fit  errors  were  in  excess  of  3000%.  Poor  fits  were  almost 
exclusively  due  to  the  type  of  error  function  used  to  train  the  network. 

As  mentioned,  training  is  based  on  minimizing  the  scaled  square  error.  In  addition  to 
accommodating  the  Sigmoid  function,  scaling  is  done  to  eliminate  bias  towards  larger 
variables  of  the  output  array.  Nonetheless,  the  un-scaled  absolute  error  of  a  given  fit  is 
proportional  the  spread  of  the  respective  training  variable.  Consequently,  the  relative 
error  of  small  values  can  be  substantial  for  data  with  a  large  spread.  Figure  52  shows  the 
effect  of  training  the  NN  on  a  single  variable  at  various  spread  ratios  in  two  training 
values.  As  can  be  seen,  the  network  tends  to  favor  the  larger  training  variable  as  the 
spread  increases.  This  can  become  especially  problematic  for  training  data  with  large 
outliers.  The  problem  can  be  overcome  through  the  use  of  an  error  function  based  on 
relative  error  but  will  require  re-formulation  of  the  training  algorithm. 


Figure  52:  Training  on  a  Progressively  Large  Spread  in  the  Training  Set  (Two  Values) 
Showing  a  Significant  Favor  Towards  the  Larger  Output 

Predictability  of  the  NN  was  established  by  calculating  a  prediction  risk  based  on  NCV. 
As  the  size  of  the  excluded  subsets  used  in  NCV  decrease  (in  relation  to  the  full  data  set), 
prediction  risk  reliability  increases.  Unfortunately,  with  limited  training  data,  the 
excluded  subset  had  to  be  relatively  large.  Tables  12  and  13  show  examples  of  prediction 
errors  found  by  removing  6  training  sets. 
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Table  12:  Predicted  11%  Values  Obtained  by  Removing  All  of  the  Listed  Solvents  from  the  Training  Set  (13  Hidden  Neurons  and 
Single  Output  Variable) 


E  (MPa) 

Density  (g/ccm) 

(10 

Syn. 

D0G/30min) 

(201 

Syn. 

)0G/30min) 

G’  (Pa) 

Solvent 


TRU 

E 


NitroEthan 

e 


175 

1 


Triethyl 

Amine 


100 

9 


1- 

Propanol 


117 

6 


Dibutyl 

Amine 


133 

8 


1- 

Pentanol 


142 

6 


n- 

ButylAmin 

e 


128 

9 


Predict 
ed 


%Err 

or 


TR 

UE 


Predicte 
d 


1881 


1022 


1022 


1384 


1353 


1284 


11.1 

9 


10.8 

5 


0.9 


|0.8 

8 


|0.9 

2 


1 0.8 
4 


%Err 

or 


TRU 

E 


1.22 


0.85 


0.80 


0.92 


0.86 


0.87 


117.2 

4 


137.2 

0 


142.8 

8 


145.6 

3 


149.5 

6 


1 43.1 
4 


Predicte 


17.19 


27.59 


43.07 


24.61 


23.70 


18.50 


Table  13:  Predicted  9%  Values  Obtained  by  Removing  All  of  the  Listed  Solvents  from  the  Training  Set  (13  Hidden  Neurons  and  Single 


E  (MPa) 

Density  (g/ccm) 

Syn.  (1000G/30min) 

Syn.  (2000G/30min) 

Yeild  (Pa) 

G’  (Pa) 

Solvent 


TRUE 


Predicted  %ErrorTRUEPredicted%ErrorTRUE 


Methanol 


1142 


NitroCH4 


2029 


NitroProp 


1621 


1196.75 


0.91 


1765.03 


1614.04 


1.27 


0.93 


0.94 


0.97 


Predicted  %ErrorTRUE 


46.7 


23.9 


29.0 


19.8 


36.4 


58.8 


Predicted  %ErrorTRUE 


31.6 


3.06 


9.49 


25.66 


9.90 


36.87 


1000 


1100 


0.12 


237.14 


316.22 


Predicted|%Error|TRUE 

■729 
7170 
3383 


Predicted  %Error 


1594.86 


4453.89 


381 .83 
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NitroEthane 

1723 

1692.96 

Triethyl 

Amine 

1025 

943.15 

1 -Propanol 

1324 

1029.03 

1.16 


0.82 


0.91 


0.97 


0.80 

0.88 


27.6 


43.5 

52.4 


49.2 


23.1 

26.7 


29.31 


13.55 

30.07 


437.44 


393.22 

0.10 


9893 


14300 


1000 


1200.02 


8155.96 

546.40 
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It  can  be  seen  that  predictions  of  density  and  elasticity  were  on  average  much  better  than 
all  other  predictions.  In  addition,  the  maximum  prediction  error  of  these  values  is 
significantly  less  than  variability  in  the  test  data  (i.e.  57%  in  Density  and  98%  in  E  for  the 
11%  data),  suggesting  that  the  network  is  not  merely  making  predictions  based  on  a  bulk 
threshold.  The  density  and  E  prediction  errors  are  also  very  similar  in  size.  This  is  likely 
due  to  the  inherent  similarity  between  the  variables  (i.e.  E  =  c2-p). 

One  possible  reason  that  density  and  elasticity  predictions  were  more  accurate  is  that 
their  values  were  measured  at  Huntsville  and  did  not  undergo  significant  degradation  due 
to  long  transportation  and  storage  times.  A  single  repeatability  study  performed  by  NJIT 
concluded  a  44%  change  in  syneresis  over  an  8-month  period.  This  amount  of  error  can 
significantly  warp  the  results  of  the  network.  Table  14  shows  an  example  of  the  effect 
training  set  error  has  on  elasticity  prediction. 


Table  14:  Comparison  of  Predictions  Made  with  Arbitrary  Error  (-30%)  Imposed  on  One 
of  Ten  Training  Values  (Fit  Error  Was  Within  0.001%  for  Each  Case).  Showing 


E  (MPa) 

E  (MPa) 

trained  with  added 

error 

Solvent 

TRUE 

Predicted 

NitroEthane 

1751 

1881 

Triethyl  Amine 

1009 

1022 

1 -Propanol 

1176 

1022 

Dibutyl  Amine 

1338 

1384 

1-Pentanol 

1426 

1353 

n-ButylAmine 

1289 

1284 

Average 


%  E rro r  Predicted  %  E rro r 
1708 


1268 


923 


1533 


1170 


1791 


20.2 


In  Table  14,  1  of  10  training  values  was  arbitrarily  reduced  by  30%  and  the  network  was 
re-trained  to  predict  the  same  set  of  values  as  the  undisturbed  case.  A  reduction  was  used 
to  insure  that  the  influence  of  variable  spread  was  not  introduced.  As  can  be  seen,  a 
relatively  modest  error  in  just  one  of  ten  training  values  can  significantly  alter  the 
networks  prediction  ability.  The  average  prediction  error  rose  from  5.1  to  20.2%.  An 
average  error  of  only  9%  was  obtained  by  leaving  the  training  set  (associated  with  the 
same  reduced  value)  out  of  training  entirely;  suggesting  that  missing  data  is  less 
detrimental  than  erroneous  data. 

As  apposed  to  merely  removing  an  entire  data  set  or  variable  from  training,  due  to  a 
single  missing  or  suspected  faulty  value  (or  even  worse  supplying  arbitrary  values  to 
satisfy  missing  data),  CFDRC  elected  to  update  the  error  function  so  that  absent  data  is 
simply  ignored.  In  this  way,  the  network  utilizes  the  most  relevant  information  available, 
which  is  critical  with  limited  training  sets.  Better  predictability,  on  data  with  inherent 
error,  might  be  achieved  in  the  future  through  larger  training  sets  involving  repeat 
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samples.  The  network  should  disregard  random  error  in  the  data  if  the  number  of  repeats 
is  sufficient  enough  to  establish  the  error  as  random  (i.e.  un-correlated). 

As  mentioned  earlier,  the  training  set  spread  has  an  unfavorable  influence  on  the 
networks  ability  to  form  accurate  fits.  Consequently,  the  upper  limit  of  predictability  will 
suffer  for  training  data  with  a  large  spread.  In  order  to  achieve  better  accuracy  without 
modifying  the  existing  training  algorithm,  a  procedure  was  established  where  the  network 
is  re-trained  on  the  fit  error.  The  new  “error  network”  (  NNerror )  will  favor  large  values 

and  accordingly  counter  balance  the  effect  of  large  spread  ratios  in  the  original  network. 
The  outputs  of  the  two  networks  are  then  added  for  each  input  array  (7)  to  obtain  mi 
overall  output  (O): 


NN(I)  +  NNerr0r(I)  =  O 

It  is  noted  that  the  network  size  will  effectively  grow  during  this  process,  but 
significantly  outperformed  (in  terms  of  speed  and  accuracy)  a  single  network  with 
comparable  dimensions.  In  addition,  the  error  network  is  not  required  to  be  the  full  size 
of  the  originally  trained  network.  Figure  53  shows  the  effect  of  this  procedure  while 
training  on  single  output  variables  of  the  11%  database.  The  procedure  effectively 
broadens  the  networks  ability  to  look  past  the  first  minimum  it  encounters  in  training 
error,  which  is  likely,  why  it  outperformed  a  single  network  of  much  larger  size.  This  is 
due  to  the  fact  that  the  first  error  minimum  reached  in  training  is  not  necessarily  optimal 
and  correlation  potential  may  still  exist  (inside  the  error)  at  this  point.  As  such,  the 
procedure  may  be  repeated  indefinitely  or  until  the  fit  is  truly  optimal. 
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20.0 


Figure  53:  Comparison  of  Average  Fit  Errors  in  a  Re-Trained  “Error  NN”  to  the 
Original  NN  Indicating  Significant  Potential  to  Increase  Overall  Predictive  Quality 

Another  issue  that  was  addressed  is  connected  to  the  concept  of  over-fitting.  If  the  error 
is  not  random  (i.e.  established  as  random  through  significant  repeat  data),  a  large  network 
will  tend  to  fit  local  correlations  in  the  error.  These  correlations  can  be  magnified  in  the 
predictions  of  new  data.  In  addition,  data  that  does  not  have  an  adequate  distribution  can 
create  similar  results.  To  address  these  problems,  SNC  was  used  to  establish  an  optimal 
network  size  (number  of  hidden  neurons).  Figure  54  shows  the  path  taken  by  SNC  in 
establishing  the  optimal  number  of  hidden  neurons  while  training  on  a  single  output 
variable.  As  can  be  seen  a  larger  network  doesn’t  necessarily  correspond  smaller 
prediction  risk  (PR9  <  PR11),  although,  it  does  correspond  to  the  smaller  square  error. 
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Figure  54:  Path  Taken  by  SNC  in  Determining  the  Optimal  Number  of  Hidden  Neurons 

(1 1%-G’  Data) 

The  benefit  of  using  SNC  goes  beyond  finding  an  optimal  number  of  hidden  neurons  to 
reduce  prediction  risk.  As  can  be  seen  in  Figure  55,  SNC  can  reduce  the  prediction  risk 
associated  with  a  given  network  size  altogether. 
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Figure  55:  Comparison  of  Networks  ( Same  Size)  Trained  with  SNC  and  Random 
Initialization.  Indicating  That  SNC  has  Potential  Beyond  Finding  an  Optimal  Number  of 

Hidden  Neurons 


Figure  55  compares  the  training  results  of  SNC  produced  networks  and  their  random  start 
counterparts.  As  can  be  seen  the  SNC  network  tends  to  converge  to  an  overall  better 
prediction  risk  than  training  a  network  from  a  random  weight  initialization.  This  takes 
place  because  as  the  SNC  network  grows;  it  accumulates  the  traits  of  the  networks  before 
it.  Starting  from  a  random  initialization  with  a  large  number  of  neurons  will  likely  result 
in  the  network  falling  on  a  minimum  error  sooner  than  in  SNC,  but  with  a  larger  value. 
Consequently,  the  networks  should  be  periodically  optimized  using  SNC  when  new 
training  data  becomes  available.  Otherwise,  the  networks  can  be  re-trained  by  initializing 
the  weights  to  those  last  obtained  by  SNC. 
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4.3.3  Python  GUI 

The  primary  objective  of  this  task  was  to  design  and  implement  a  graphical  user  interface 
(GUI)  for  the  NN  algorithms  and  gel  database  files.  The  GUI  facilitates  NN  training, 
prediction,  and  database  organization  through  a  windows-like  environment. 
Consequently,  user  understanding  of  the  underlying  principles  and  details  is  not  required. 
In  addition,  the  GUI  can  be  loaded  on  any  platform. 

4.3.3.1  Graphical  User  Interface  Application  Development 

The  Python  programming  language  (http://www.python.org/)  was  selected  for  the 
implementation  of  the  GUI.  Python  is  an  open  source,  platform  independent 
programming  language  with  minimal  syntax  and  semantics,  facilitating  rapid 
implementation.  The  Python  standard  library  is  comprehensive,  providing  support  for 
web  and  GUI  application  development.  The  Python  standard  library  was  supplemented 
with  the  wxPython  library,  an  open  source,  platform  independent  GUI  library  developed 
from  the  wxWidgets  library  and  the  Python  programming  language.  The  wxPython 
library  enables  the  development  of  applications  that  utilize  the  controls  and  utilities  of  the 
native  platform,  in  this  case  Microsoft  Windows.  Extensible  markup  language  (XML) 
read  and  write  functions  were  implemented  utilizing  the  open  source  pyXML  library. 

4.3.3.2  Development  of  a  Chemical  Compound  Library 

A  chemical  compound  library  was  developed  utilizing  XML.  A  primary  advantage  to 
using  XML  is  that  information  may  be  described  in  a  hierarchical  manner,  facilitating  its 
proper  annotation.  The  chemical  compound  library  contains  physical  properties 
information  for  three  classes  of  compounds:  gellants,  surfactants  and  liquids.  The  custom 
XML  schema  of  the  chemical  compound  library  is  illustrated  in  figure  56.  The  chemical 
compound  library  was  populated  with  physical  property  information  from  Tables  7  and  8, 
which  includes  compound  name,  CAS  number,  molar  mass,  density,  melting  point, 
boiling  point,  standard  enthalpy  change  of  formation,  viscosity,  surface  tension,  dielectric 
constant,  vapor  pressure,  accentric  factor  and  dipole  moment.  In  addition,  physical 
properties  associated  with  gel  compounds  are  included  in  the  chemical  compound  library. 
These  physical  properties  include  density,  syneresis,  storage  modulus,  bulk  modulus  and 
50-caliber  test.  At  the  moment,  the  physical  properties  included  in  the  chemical 
compound  library  enable  the  calculation  of: 

□  Storage  Modulus 

□  Yield  Stress 

□  Syneresis 

o  2000G/90min 
o  100000G/30min 

□  Density 

□  Modulus  of  Elasticity 

Calculation  of  other  properties  can  be  made  possible  through  a  slight  restructuring  of  the 
GUI. 
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<?xml  version="1.0"  encoding="UTF-8"  ?> 

-  <compoundl_ibrary> 

-  <listOfGellants> 

-  <gellant  name="Sodium  silicate"  CASNumber="6834-92-0"> 

<molarMass  unit="g/mol"  value="122.06"  /> 

<density  unit="g/cm3"  value="2.4"  /> 

</gellant> 

</listOfGellants> 

-  <listOfSurfactants> 

-  <surfactant  name="Sodium  dodecyl  sulfate"  CASNumber="151-21-3"> 

<molarMass  unit="g/mol"  value="288.38"  /> 

<density  unit="g/cm3"  value="1.01"  /> 

</surfactant> 

</listOfSurfactants> 

-  <listOfLiquids> 

-  <liquid  name=' iso-Pentane"  CASNumber="78-78-4"> 

<molarMass  unit="g/mol"  value="72.15"  /> 

<density  unit="g/cm3"  value="0.6"  /> 

<meltingPoint  unit="K"  value="-160.0"  /> 

<boilingPoint  unit="K"  value="27.88"  /> 

<triplePoint  unit="K"  value="113.37"  /> 

<criticalTemperature  unit="K"  value="460.35"  /> 
ccriticalPressure  unit="kPa"  value="3395.7"  /> 

<standardEnthalpyChangeOfFormation  unit="kJ/mol"  value="-179.0"  /> 
<standardMolarEntropy  unit="J/molK"  value="260.7"  /> 

<hedlCdpduily  uriil="J/rnulK"  vdlue="164.85"  /> 

<viscosity  unit="mPas"  value="0.214"  /> 

<surfaceTension  unit="mN/m"  value="14.45"  /> 

<dielectricConstant  value="  1.845"  /> 

<thermalConductivity  unit="w/mK"  value="0.111"  /> 

<vaporPressure  unit="C"  value="79.3"  /> 

<acentricFactor  value="0.2296"  /> 

<dipoleMoment  unit="D"  value="0.1"  /> 

</liquid> 

-  diquid  name="Hexene"  CASNumber="592-41-6"> 

<molarMass  unit="g/mol"  value="84.1608"  /> 

<density  unit="g/cm3"  value="0.6731"  /> 

<meltingPoint  unit="K"  value="-140.0"  /> 
cboilingPoint  unit="K"  value="63.48"  /> 

<standardEnthalpyChangeOfFormation  unit="kJ/mol"  value="-74.2"  /> 

<viscosity  unit-'mPas"  value="0.252"  /> 

<surfaceTension  unit="mN/m"  value="17.9"  /> 

<dielectricConstant  value="2.077"  /> 
cvaporPressure  unit="C"  value="24.8"  /> 

<acentricFactor  value="0.281"  /> 

<dipoleMoment  unit="D"  value="0"  /> 

</liquid> 

</listOfLiquids> 

_ </compoundLibrary> _ 

Figure  56:  The  custom  XML  schema  of  the  chemical  compound  library 


4.3.3.3  Design  and  Implementation  of  Chemical  Compound  Library  Management 
Page 

The  chemical  compound  library  management  page  was  designed  and  implemented  to 
enable  a  user  to  view,  create,  edit  and  delete  gellant,  surfactant  or  liquid  compounds  in 
the  chemical  compound  libraiy. 

A  user  may  navigate  to  the  ‘compound  library’  page  using  the  page  tabs  at  the  top  of  the 
application  window  (Figure  57).  The  default  ‘compound  class’  is  liquid.  A  user  may 
select  the  gellant,  surfactant  or  gel  radio  buttons  in  the  compound  type  panel  in  order  to 
view  the  gellant,  surfactant,  or  gel  compounds  and  associated  physical  properties.  The 
compound  names  associated  with  the  selected  compound  class  existing  within  the  library 
are  displayed  in  the  list  box  within  the  ‘compound  name’  panel.  A  user  may  select  a 
compound  name  to  view  the  physical  properties  associated  with  the  compound  name 
existing  within  the  library.  A  user  may  create  a  new  compound  by  selecting  the  desired 
compound  class,  clicking  the  ‘new’  button  and  entering  the  associated  physical  property 
information  in  the  text  boxes  within  the  ‘compound  properties’  panel.  A  user  may  also 
edit  an  existing  compound  by  selecting  the  desired  compound  class,  selecting  the  desired 
compound  name  and  editing  the  associated  physical  property  information.  Finally,  a  user 
may  delete  an  existing  compound  by  selecting  the  desired  compound  class,  selecting  the 
desired  compound  name,  and  clicking  the  ‘delete’  button.  A  warning  dialog  box  will  be 
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displayed  and  the  user  may  click  the  ‘yes’  button  to  delete  the  compound  or  the  ‘no’ 
button  to  return  to  the  compound  library  page.  A  user  may  click  the  ‘clear’  button  to 
remove  compound  name  selection  and  associated  physical  properties. 

The  chemical  compound  library  management  page  was  designed  and  implemented  to 
include  a  minimal  and  intuitive  interface  in  order  to  reduce  the  burden  on  the  user.  In 
addition,  the  user  is  aided  by  error  messages  displayed  for  nonnumeric  physical  property 
entries.  Limiting  user  error  in  the  management  of  the  chemical  compound  library  will 
assist  in  limiting  error  in  calculation. 


Figure  57:  Screenshot  of  chemical  compound  library  management  page 

4.3.3.4  Design  and  Implementation  of  Calculation  Page 

The  calculation  page  was  designed  and  implemented  to  enable  a  user  to  calculate  the 
physical  properties  of  a  simulated  propellant.  The  ‘calculate’  page  is  the  default  page  of 
the  appli cation  window  (Figure  58),  however,  a  user  may  navigate  to  the  calculate  page 
using  the  page  tabs  at  the  top  of  the  application  window.  In  order  to  perform  a 
calculation,  a  user  must  select  a  compound  name  from  the  ‘compound  name’  panel.  The 
user  must  then  select  the  training  library  file  used  in  the  calculation  by  entering  the  file 
location  in  the  text  box  within  the  ‘training  library’  panel  or  by  clicking  the  browse 
button  and  selecting  the  file  from  the  ‘Choose  a  file’  dialog  box.  In  order  to  perform  the 
calculation  with  the  selected  compound  name  and  training  library  file,  the  user  must  click 
the  ‘calculate’  button.  The  user  may  clear  the  entered  information  by  clicking  the  ‘clear’ 
button.  Results  of  the  calculation  are  displayed  in  the  ‘results’  panel. 
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4.3 .3.5  Design  and  Implementation  of  Training  Page 

The  training  page  was  designed  and  implemented  to  enable  a  user  to  select  multiple  gels, 
available  from  the  chemical  compound  library,  to  view  their  associated  physical 
properties  and  to  calculate  a  new  training  set,  which  may  be  saved  by  the  user. 

A  user  may  navigate  to  the  ‘training’  page  using  the  page  tabs  at  the  top  of  the 
application  window  (Figure  58).  A  user  may  then  select  one  or  more  gel  names  from  the 
list  box  within  the  ‘gel  name’  panel  and  click  the  ‘add’  button  to  append  the  gel  name  and 
its  associated  physical  properties  to  the  training  list.  A  user  may  select  a  gel  name  in  the 
training  list  and  click  the  ‘remove’  button  to  remove  the  gel  name  and  its  associated 
physical  properties  from  the  training  list.  A  user  may  remove  all  gel  names  and 
associated  physical  properties  by  clicking  the  ‘clear’  button.  In  order  to  create  a  new 
training  set,  a  user  may  add  the  desired  gel  names  and  associated  physical  properties  to 
the  training  list  and  click  the  ‘train’  button.  The  results  of  the  training  may  be  saved  to  a 
desired  location  by  clicking  the  ‘save’  button. 


File  Help 

Calculate  Train  |  Compound  Library  | 
Training  Properties 
Propellant  Name 

Benzere  |  9%  Silica 
Cyclohexane  |  9%  Silica 
Cyclohexene  j  9%  Silica 
Decane  |  9%  Silica 
Decanol  |  9%  Silica 


^]n]xj 


File  Help 

Calculate  |  Train  |  Compound  Library  | 
Compound  Name 


-=JnJ2<] 


Acetone 

Cyclohexane 

Cyclohexene 

d 

Training  Library 
File  |  | 


Figure  58:  SOGeF  GUI  Application  Windows:  (Left)  Screenshot  of  calculate  page 

(Right)  Screenshot  of  training  page 
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4.4  Dynamic  Brownian  /  Stokesian  Gelation  Simulation 


The  modeling  techniques  of  Brownian  and  Stokesian  dynamics  were  established  for 
strong  attractive  forces  in  Phase  I.  The  focus  of  the  Phase  II  effort  was  to  develop 
parameterized  forces,  which  include:  van  der  Waals,  solvation,  depletion,  and  repulsive 
charge.  Unfortunately,  these  forces  depend  on  the  type  of  solvent  used.  To  date,  only 
water  solvents  have  been  modeled. 

In  the  Brownian  dynamics  approach,  the  silica  particles  are  assumed  to  be  much  larger 
than  the  solvent  particles.  This  approach  is  justified,  in  our  case,  by  the  fact  that  the  silica 
particles  are  estimated  to  be  at  least  40,000  times  more  massive  than  the  solvent 
molecules.  Consequently,  the  solvent  was  modeled  as  a  stationary  continuum  and  only 
the  motion  of  the  silica  particles  were  tracked.  The  particles  move  by  the  influence  of 
inter-particle  forces  as  well  as  the  drag  exerted  by  the  solvent.  In  addition,  random 
thermal  motion  is  specified  through  Brownian  modeling. 

The  inter-particle  force  potential  ( U )  was  modeled  with  the  following  equation: 

=  4 A  ■  [(2  /  r )12  -  (2  /  r )6  +  exp (-r  /  2)  /  r ] 

KT 

Where  r  is  the  non-dimensional  distance  (with  respect  to  particle  radius)  between  particle 
centers  (i.e.  r  =  2  implies  contact).  Figure  59  shows  a  potential  barrier  occurring  just 
before  contact  followed  by  a  closer  stable  “bonding”  well. 


Figure  59:  Form  of  the  inter-particle  force  typical  of  silica  in  aqueous  environments.  The 
repulsive  portion  with  local  maximum  away  from  contact  (r=2)  is  due  to  charged  surface 
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Figures  60  and  6 1  shown  comparisons  of  the  Stokesian  and  Brownian  models  for  the 
prediction  of  mean  square  displacement  and  number  of  particle  bonds  respectively. 
Notice  that  bonds  can  be  broken  (Figure  61)  and  that  the  Stokesian  model  shows  a  slower 
gelation  process  (due  to  a  larger  hydrodynamic  force). 


Figure  60:  Mean  square  displacement  (in  units  of  particle  radius  a)  as  a  function  of 
dimensionless  time,  t*  =  Dt/a2  with  D  the  bore-particle  diffusion  constant  ( Stokes - 
Einstein)  for  simulations  by  Brownian  and  Stokesian  Dynamics  (BD,  SD)  of  identical 
initial  configurations ;  the  simulations  differ  in  BD  lacking  hydrodynamic  interactions 
which  are  rigorously  accounted  in  SD.  These  are  conditions  yielding  strong  gels 
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Figure  61:  Number  of  bonds,  defined  as  particle  pairs  having  centers  within  2.3  radii,  as 
a  function  of  dimensionless  time  with  124  particles  in  the  unit  cell 

One  flaw  in  the  model  is  that  the  predicted  aggregates  form  densely  packed  clusters  (see 
Figure  62:  Left),  as  apposed  to  a  space-filling  network  typical  of  gels.  The  flaw  is 
quantified  by  the  pair  distribution  shown  in  Figure  62  (Right).  Notice  the  high  number  of 
particles  inside  the  “bonding”  well.  This  tendency  is  believed  to  occur  as  a  result  of 
inadequate  contact  “frictional”  modeling,  which  would  prevent  the  fractal  structures  from 
collapsing. 


Figure  62:  (Left)  Snapshot  of  particle  configuration  for  an  initially  uniformly  distributed 
system,  with  1000  particles  in  the  unit  cell,  at  solids  fraction  of  (j)  =  0.05  and  the  force 
parameter  (see  Fig  3)  of  A  =20.  (Right)  The  pair  distribution  function  of  a  gelling 
system,  showing  large  value  at  contact  and  local  crystalline  ordering. 
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The  simulation  tools,  developed  in  Phase  I  mid  II,  have  two  major  drawbacks.  First,  the 
Accelerated  Stokesian  Dynamics  model  (ASD  as  described  in  10)  limits  system  size  and 
requires  excessive  run-time  to  achieve  gelation,  due  to  rigorous  hydrodynamic 
calculations.  Secondly,  the  Brownian  Dynamics  (BD)  model,  while  allowing  much 
larger  system  sizes  and  reduced  computation  time,  does  not  account  for  realistic 
hydrodynamic  interaction.  A  compromise  can  be  established  between  the  two  models,  as 
in  (11),  where  only  hydrodynamic  interactions  of  close  pair  particles  are  modeled.  Far 
pair  hydrodynamic  interactions  are  captured  through  a  “mean  field”  approach,  which 
involves  a  systematic  account  of  the  added  dissipation.  This  accelerates  the  overall 
method  by  more  than  mi  order  of  magnitude,  with  little  loss  in  physical  accuracy. 
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5. 


SUMMARY 


The  main  goal  of  this  STTR  Phase  II  project  was  to  develop  a  software  tool  to  assist  in 
the  synthesis  and  optimization  of  gelled  propellants.  CFDRC  successfully  accomplished 
this  goal  with  the  creation  of  SOGeF.  SOGeF  is  a  neural  network-based  QSPR  gel 
prediction  algorithm  wrapped  in  a  convenient  and  robust  GUI.  The  main  advantage  of 
SOGeF  over  physical-based  models  is  that  input  parameters  are  easily  identified  at  the 
macroscopic  level.  In  addition,  SOGeF  predictions  are  explicit  and  do  not  require 
lengthy  gelation  simulations. 

Validation  and  training  of  SOGeF  began  with  the  creation  of  an  extensive  gel  database. 
In  Phase  II,  over  39  different  organic-based  gels  were  prepared,  whose  various  physical 
properties  were  measured  by  CFDRC,  NJIT,  mid  CCNY.  Comparisons  of  predictions 
made  by  SOGeF  to  removed  elements  in  the  training  database  were  promising  and  it  is 
expected  that  accuracy  will  only  get  better  as  the  number  of  training  sets  grows. 
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6. 


CONCLUSIONS 


The  overall  conclusion  from  the  Phase  II  study  is  that  SOGeF  is  an  extremly  powerful 
tool  for  the  synthesis  and  optimization  of  gelled  propellants.  A  key  advantage  of  SOGeF 
is  that  model  predictions  are  explicit  and  do  not  require  lengthy  gelation  simulations.  In 
addition,  SOGeF  input  parameters  are  easily  identified,  as  apposed  to  the  challenging 
specitication  physical  identifiers  in  gelation  models.  Also,  as  the  size  of  the  training 
database  grows,  predictions  from  SOGeF  will  only  become  more  reliable,  wheras  the 
accuracy  of  physical  models  is  relatively  fixed. 

The  main  conclusions  of  the  Phase  II  work  can  be  summarized  as  follows: 

1.  Gel  preparation 

a.  It  is  necessary  to  mix  gels  thoroughly  mid  relatively  quick.  Once  a  gel  is 
established,  it  is  impossible  to  increase  its  silica  content. 

b.  Preparation  of  alcohol -based  gels  requires  the  use  of  surfactants  in  order  to 
reduce  interfering  hydrogen  bonds  between  the  solvent  and  silica  network. 

c.  Strong  gels  were  formed  from  alkanes,  alkenes,  mid  cyclic  alkanes,  while 
weak  gels  were  formed  from  the  alcohols  and  acetone  solvents. 

d.  The  influence  of  operator  variability  in  preparing  the  gels  had  very  little 
impact  on  their  consistency. 

e.  It  was  necessary  to  increase  the  silica  content  of  the  organic  gels  from  9  to 
1 1%  in  order  to  improve  the  50-caliber  test  pass  rate. 

2.  Measurement  observations 

a.  Silica  content  has  the  most  influence  on  the  gel  properties  of  a  specific 
organic  group. 

b.  Storage  time  has  a  substantial  impact  on  the  properties  of  organic-based 
gels.  This  is  believed  to  be  a  result  of  solvent  evaporation. 

c.  The  syneresis  of  organic-based  gels  shows  a  strong  dependence  on  the 
magnitude  of  acceleration  exposure,  with  a  minimal  dependence  on  the 
duration  of  exposure.  Once  a  gel  is  exposed  to  a  given  acceleration,  it 
immediately  loses  a  fixed  amount  of  solvent. 

d.  The  presence  of  an  -OH  functionality  in  a  solvent  results  in  poor  gelation. 

e.  DLS  and  rotational  rheometry  form  comparable  estimates  of  gelation  time. 

f.  Salt  was  critical  to  the  formation  of  aqueous  colloidal  silica  gels,  as  it 
destabilizes  the  repulsive  silica  charges.  The  effect  becomes  saturated  at 
large  salt  levels,  where  the  structural  growth  is  limited  solely  by  diffusion. 

3.  Neural  Network  functionality 

a.  The  chosen  descriptor  set  forms  well-defined  PDI  patterns  in  the  alcohol- 
based  gels.  This  was  mostly  a  result  of  the  inherent  ability  of  accentric 
factor  to  distinguish  organic  groups. 

b.  In  general,  error  in  the  training  fit  forms  an  upper  bound  on  predictability. 

c.  Spread  in  the  training  data  directly  scales  into  the  errors  seen  in  the 
training  fit.  Consequently,  relative  error  will  be  biased  towards  larger 
values  in  the  training  set. 
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d.  Re-training  the  network  on  its  fit  error  can  counteract  issues  arising  from 
training  variable  spread  and  local  error  minimums.  The  size  of  the  re¬ 
trained  network  can  be  much  smaller  than  the  original  network. 

e.  Errors  in  training  data  can  significantly  warp  the  networks  output. 
Consequently,  it  is  crucial  that  missing  data  not  be  replaced  by  arbitrary 
values  (the  network  was  designed  to  disregard  missing  data). 

f.  SNC  has  value  not  only  in  selecting  an  appropriate  network  size,  but  also 
in  reducing  the  prediction  risk  associated  with  that  size.  It  is  therefore 
advisable  to  periodically  re-train  the  network  using  SNC  as  new  training 
data  becomes  available. 
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RECOMMENDATIONS 


7  . 


□  Gel  Database 

o  Collect  more  data  to  broaden  the  networks  prediction  ability 

o  Collect  repeat  samples  to  establish  random  error  in  network  training 

□  Neural  Network  ( QSPR ) 

o  Further  optimize  descriptor  variable  selection,  as  an  automated  process 
(Input  Variable  Pruning) 

o  Explore  the  use  of  different  error  functions  (relative  error)  to  eliminate 
variable  spread  issues 

o  Detect  and  optimize  the  networks  ability  to  establish  random  error  in  the 
training  data 

o  Further  explore  network  error  re-training  to  obtain  optimal  models 

o  As  new  training  data  is  acquired,  the  network  should  be  periodically  re¬ 
trained  using  SNC 

□  Gelation  Modeling 

o  Determine  force  fields  for  various  organic  solvent  types 

o  Adopt  a  compromised  model,  as  in  (11),  to  decrease  run-time  and  increase 
system  size 

o  Include  particle-particle  friction  so  that  space  filling  networks  can  be 
modeled  (apposed  to  the  current  collapsed  networks) 


73 


8731/Final  Draft 


9. 
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Appendix  A  (Solvent  Descriptor  Values) 


lOOmsc 


(g/ccm) 

25°C 

(mPa*sec) 

25°C 

(mN/m) 

25°C 

(20°  C) 

kPa 

(25°C) 

(Debye) 

(kJ/mol) 

298K 

(g/1 00m  L 
) 

(°C) 

(°C) 

Comp/Param 

Rho 

Vise 

SurfTens 

DielCon 

st 

VapPress 

AccFact 

Dipole 

M 

Std  HF 

WaterSol  MeltP 

BoilP 

iso-Pentane 

0.6 

0.214 

14.45 

1.845 

79.3 

0.2296 

0.1 

-179 

0 

-160 

27.88 

Heptane 

0.684 

0.387 

19.66 

1 .9209 

6.09 

0.349 

0 

-224.4 

0 

-91 

98.4 

Octane 

0.703 

0.508 

21.14 

1.948 

1.86 

0.393 

0.07 

-250.3 

0 

-57 

125.67 

iso-Octane 

0.688 

0.473 

18.77 

1.94 

6.6 

0.303 

0 

-259 

0 

-107 

99.22 

CycloHexane 

0.7785 

0.894 

24.65 

2.0243 

13 

0.211 

0 

-156.4 

0 

7 

80.73 

Hexene 

0.6731 

0.252 

17.9 

2.077 

24.8 

0.281 

0 

-74.2 

0.005 

-140 

63.48 

Benzene 

0.8765 

0.604 

28.22 

2.2825 

12.7 

0.209 

0 

49.1 

0.179 

6 

80.09 

Diethyl  Amine 

0.7056 

0.319 

19.85 

3.68 

30.1 

0.291 

0.92 

-103.7 

81 

-50 

55.5 

n-PropylAmine 

0.7173 

0.376 

21.75 

5.08 

42.1 

0.303 

1.17 

-101.5 

100 

-83 

47.22 

Cyclohexene 

0.811 

0.625 

26.17 

2.2176 

11.8 

0.2312 

0.332 

-38.5 

0 

-103 

82.98 

n-Pentane 

0.626 

0.224 

15.49 

1.8371 

68.3 

0.251 

0 

-173.5 

0.01 

-130 

36.06 

Hexane 

0.6548 

0.294 

17.59 

1 .8865 

20.2 

0.299 

0.05 

-198.7 

0.0013 

-95 

68.73 

Decane 

0.73 

0.838 

23.37 

1 .9853 

0.2 

0.488 

0.07 

-301 

0 

-30 

174.15 

MEK 

0.8054 

0.405 

23.97 

18.56 

12.6 

0.32 

2.78 

-273.3 

29 

-87 

79.59 

Acetone 

0.7899 

0.306 

23.46 

21.01 

30.8 

0.306 

2.88 

-248.4 

100 

-95 

56.05 

MIBK 

0.7978 

0.545 

23.9 

13.11 

2.64 

0.37 

4.2 

-305.3 

1.91 

-84 

116.5 

Decanol 

0.8297 

10.91 

28.51 

7.93 

0.009 

0.624 

1.8 

-478.1 

0.37 

7 

231.1 

Methanol 

0.7914 

0.544 

22.07 

33 

16.9 

0.565 

1.7 

-239.2 

100 

-97 

64.6 

NitroCH4 

1.1371 

0.63 

36.53 

37.27 

4.79 

0.31 

3.46 

-112.6 

10 

-28 

101.19 

NitroProp 

0.9961 

0.798 

29.29 

24.7 

1.36 

0.376 

3.66 

-167.2 

1.4 

-108 

131.1 

NitroEthane 

1 .0448 

0.688 

32.13 

29.11 

2.79 

0.341 

3.23 

-143.9 

4.6 

-89 

114 

Triethyl  Amine 

0.7275 

0.347 

20.22 

2.418 

7.7 

0.32 

0.61 

-127.7 

17 

-115 

89 

1 -Propanol 

0.8035 

1.945 

23.32 

20.8 

0.154 

0.628 

1.55 

-302.6 

100 

-126 

97.2 

surfactant 


surfactant 
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Dibutyl  Amine 

0.767 

0.918 

24.12 

2.765 

0.34 

1 

-206 

3.5 

-62.15 

159.6 

1-Pentanol 

0.8144 

3.5 

25.6 

13.9 

0.293 

0.621 

1.712 

-295.63 

2.1956 

-78.2 

137.983 

n-ButylAmine 

0.74 

0.5 

16.87 

5.4 

12.23 

0.285 

1.361 

-94.98 

100 

-49 

77 
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Appendix  B  (Additional  Syneresis  Data) 


Solvent/Liquid 


Serial  No. 


Composition 


5min  90min 


n-Pentane 


Iso-Pentane 


Hexane 


Heptane 


Octane 


Iso-Octane 


Decane 


082307-01 


082307-02 


082207-01 


9%  Cabosi 


9%  Cabosi 


9%  Cabosi 


081407-01 

9%  Cabosil 

0.06 

0.135 

082807-01 

9%  Cabosil 

0.06 

0.134 

090607-01 

9%  Cabosil 

1.73 

1.54 

082907-01 

9%  Cabosil 

0.067 

0.73 

Solvent/Liquid 


Serial  No. 


Composition 


lumin  30min  90min 


1 

n-Pentane 

082307-01 

9%  Cabosil 

0.4 

2 

Iso-Pentane 

082307-02 

9%  Cabosil 

0.61 

Hexane 


Heptane 


Octane 


Iso-Octane 


Decane 


082207-01 


081407-01 


082807-01 


090607-01 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


082907-01 

9%  Cabosil 

032508-02 

9%  Cabosil,  no  Duomeen 

9 

|  Nitroethane  | 

022508-01 

9%  Cabosil 

10 

022708-01 

9%  Cabosil 

11 

Benzene 

021808-01 

9%  Cabosil 

12 

Acetone 

031108-01 

9%  Cabosil 

Nitromethane 


Methanol 


Cyclohexane 


1  -  Hexene 


1  -Oct  an  o  I 


1  -Butanol 


Methyl -Ethyl  Ketone 


Cyclohexene 


1-Hexanol 


Decanol 


022808-01  9%  Cabosil 


031208-01  9%  Cabosil,  0.9%  Surfactant 


021108-01  9%  Cabosil 


021908-01  9%  Cabosil 


032008-01  9%  Cabosil,  0.9%  Duomeen 


032708-01  9%  Cabosil,  0.9%  Duomeen 


032608-02  9%  Cabosil 


031308-01 

9%  Cabosil,  0.9%  Duomeen 

021308-01 

9%  Cabosil 

040208-02 

9%  Cabosil,  0.9%  Duomeen 

042208-01 

9%  Cabosil 

042408-01 

9%  Cabosil 

n-Propyl  Amine 


Hexane 


Diethyl  Amine 


Methyl -Ethyl  Ketone 


MIBK 


Nitromethane 


Decane 


Acetone 


MIBK 


Nitropropane 


Nitropropane 


Methanol 


Pentane 


042908-01 


040308-01 


043008-01 


051408-01 


052108-01 


060408-01 


050508-01 


051508-01 


052008-01 


060508-01 


060908-01 


060308-01 


050708-01 


9%  Cabosil 


11%  Cabosil 


9%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


9%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


35.64257 

38.05896 

39.10955 

4.633094 

7.477567 

7.677275 

6.028267 

7.954173 

9.494894 

5.826636 

3.969884 

7.205179 

28.56252 

29.54695 

34.17266 

1.945412 

4.154495 

3.063759 

28.5155 

29.89409 

31.6416 

8.014049 

9.330223 

10.18804 

8.090379 

6.284737 

8.571429 

6.430077 

34.1629 

32.65993 

31.96605 

35.98039 

38.34163 

39.16376 

41.15409 

40.53167 

23.31774 

34.13249 

37.04553 

1 .739467 

2.091768 

2.217529 

41.11963 


12.96417 


25.34562 


41.04596 


7.177211 


25.54897 


33.07692 


29.93241 


0.03308 


5.804579 


27.21722 


33.95739 


2.570778 


0.032478 


32.69105 


0.119713 
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s/n 


Serial  No. 


5000G 


n-Pentane 


Iso-Pentane 


Hexane 


Octane 


Iso-Octane 


Decane 


082307-01 


082307-02 


082207-01 


081407-01 


082807-01 


090607-01 


082907-01 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


1.803607 


7.138001 


9.58445 


9.986684 


16.70007 


19.1689 


18.67835 


0.475543 


6.875834 


9.716599 


10.13423 


16.19813 


21.45722 


20.52314 


2.93725 


8.198925 


9.479306 


9.152086 


16.34031 


20.09346 


17.38255 


1.615074 


10.35404 


9.659864 


11.13333 


15.98934 


20.0133 


19.09212 


3.686327 


8.916324 


12.64446 


11.32708 


18.27815 


21.9661 


21 .30376 


3.459747 


10.73826 


12.20499 


11.16342 


18.28533 


21.70963 


20.63385 


2.374491 


10.79812 


13.34675 


12.56684 


16.84492 


20.73579 


20.10688 


Solvent/Liquid 


Composition 


5min  10min  15min  20min  25min  30min 


n-Pentane 


Iso-Pentane 


Hexane 


Heptane 


Octane 


Iso-Octane 


Decane 


1  -  Propanol 


Nitroethane 


Nitropropane 


Benzene 


Acetone 


Nitromethane 


Methanol 


Cyclohexane 


1  -  Hexene 


1-Octanol 


1 -Butanol 


Methyl-Ethyl  Ketone 


1 -Propanol 


Cyclohexene 


1-Hexanol 


Decanol 


Hexane 


Diethyl  Amine 


Methyl-Ethyl  Ketone 


MIBK 


Nitromethane 


Decane 


Acetone 


MIBK 


Nitropropane 


Nitropropane 


Methanol 


Pentane 


1 -Decanol 


082307-01 


082307-02 


082207-01 


081407-01 


082807-01 


090607-01 


082907-01 


032508-02 


022508-01 


022708-01 


021 808-01 


031108-01 


022808-01 


031 208-01 


021108-01 


021 908-01 


032008-01 


032708-01 


032608-02 


031308-01 


021308-01 


040208-02 


042208-01 


042408-01 


042908-01 


040308-01 


043008-01 


051408-01 


052108-01 


060408-01 


050508-01 


051508-01 


052008-01 


060508-01 


060908-01 


060308-01 


050708-01 


060208-01 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil,  no  Duomeen 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil,  0.9%  Surfactant 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil,  0.9%  Duomeen 


9%  Cabosil,  0.9%  Duomeen 


9%  Cabosil 


9%  Cabosil,  0.9%  Duomeen 


9%  Cabosil 


9%  Cabosil,  0.9%  Duomeen 


9%  Cabosil 


9%  Cabosil 


9%  Cabosil 


11%  Cabosil 


9%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


9%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


11%  Cabosil 


5.970149 


9.538462 


13.04348 


13.23529 


19.69697 


23.8806 


22.64706 


8.296943 


12.85714 


13.63636 


4.477612 


7.575758 


18.84058 


17.46479  16.41791 


20.58824  18.46154 


24 


23.05556 


27.53623 


20.89552 


6.841339 


17.3913 


18.08346 


23.52941 


23.97661 


26.47059 


24.63768 


6.153846 


11.11111 


17.14286 


17.3913 


21.6954 


25.75758 


23.52941 


52.40642 


27.63701 


28.98507 


27.48188 


43.46532 


23.91559 


46.726 


30.58036 


31.98847 


39.55479 


47.72539 


50.50882 


47.44452 


22.769 


50.2928 


24.44366 


43.50243 


52.48933 


29.64061 


43.15899 


43.17223 


41.25989 


14.58333 


22.95835 


36.67322 


63.18153 


20.25565 


17.2391 


42.12466 


19.85637 


7.23 


79 


8731/Final  Draft 


s/n 

Solvent/Liquid 

Serial  No. 

Composition 

30000G 

5min 

90min 

1 

n-Pentane 

082307-01 

9%  Cabosil 

30.04 

29.09 

2 

Iso-Pentane 

082307-02 

9%  Cabosil 

35.4 

38.51 

3 

Hexane 

082207-01 

9%  Cabosil 

37.82 

40 

4 

Heptane 

081407-01 

9%  Cabosil 

35.66 

39.98 

5 

Octane 

082807-01 

9%  Cabosil 

38.23 

42.64 

6 

Iso-Octane 

090607-01 

9%  Cabosil 

41.58 

46.68 

7 

Decane 

082907-01 

9%  Cabosil 

36.94 

42.75 
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#include  <stdio.h> 
#include  <stdlib.h> 
#include  <math.h> 


Appendix  C  (NN  Code  Listing) 


int  Activation=l; 
int  Run=l; 

(SNC)  1 -Train  from  wo  (train  with  H  from  file) 

int  1=11; 

int  0=1; 

int  N=16; 

int  H=12; 

neurons 

char  wo_dat[50]="wo.dat",  vo_dat[50]="vo.dat"; 


//O-sigmoid  l-(x2->x) 

//O- Train  from  0 

2-Output 

//Dimension  of  input 
//Dimension  of  output 
//Number  of  training  sets 
//Maximum  #  of  hidden 


#define  maxi  1 1  //maxDimension  of 

input  (need  dynamic  memory) 

#define  maxO  6  //maxDimension  of 

output  (need  dynamic  memory) 

#define  maxN  50 
(need  dynamic  memory) 

#define  maxH  23 

hidden  neurons  (need  dynamic  memory) 

int  Itemp,  Otemp,  Ntemp,  Htemp;  //Used  to  train  subsets 

float  Input[maxN][maxI]; 
float  Output[maxN][maxO]; 

float  ScaledInput[maxN]  [maxi];  //Must  scale  unless  using  % 

error 

float  ScaledOutput[maxN][maxO]; 

float  maxIn[maxI],minIn[maxI],maxOut[maxO],minOut[maxO]; 

float  w[maxI+l][maxH],v[maxH+l][maxO];  //Weights 

including  thresholds 

float  wtemp  [maxl+ 1  ]  [maxH] ,  vtemp  [maxH+ 1  ]  [maxO] ; 

float  w_optimal[maxI+l]  [maxH], v_optimal[maxH+l]  [maxO]; 

float  Input_temp[maxN][maxI],  Output_temp[maxN][maxO]; 

int  w_test[maxI+l][maxH],v_test[maxH+l][maxO];  //set  0  to  hold 

constant 

float  mu=0.02;  //Initial  step  size 


//maxNumber  of  training  sets 
//Maximum  #  of 


void  LoadData(){ 

FILE  *LiquidFile=fopen("LiquidData.txt","r"); 
FILE  *GelFile=fopen("GelData.txt","r"); 
int  i,j,read; 
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fscanf(LiquidF  ile,  "%i"  ,&Run); 
fscanf(LiquidF  ile,  "%i"  ,& Activation); 
fscanf(LiquidF  ile,  "%i"  ,&I); 
fscanf(LiquidF  ile,  "%i"  ,&N); 
fscanf(LiquidF  ile,  "%i"  ,&H); 
fscanf(LiquidFile,"%s",&wo_dat[0]); 
fscanf(LiquidFile,"%s",&vo_dat[0]); 

fscanf(GelFile,"%i",&0); 

fscanf(GelFile,"%i",&read); 

if(read!=N){printf("Incompatible  liquid  and  gel  files\n");} 


for(i=0;i<N;i++){ 

for(j =0  ;j  <I;j  ++)  { 

fscanf(LiquidFile,"%f',&Input[i][j]); 

if(i==0){maxIn[j]=Input[i][j]*l.l;minIn|j]=Input[i]|j]*0.9;}//Larger  range  for 
extrapolation  under  sigmoid  activation 
else{ 

if(maxln[j  ]<Input[i]  [j]  *  1 .  l){maxln[j]=lnput[i]  [j]*  1.1;} 
if(minln[j]>lnput[i][j]*0.9){minln[j]=lnput[i][j]*0.9;} 

} 

} 

for(j =0  ;j  <0;j  ++)  { 

fscanf(GelFile,"%f',&Output[i][j]); 

if(i==0)  {maxOut[j  ]=Output[i]  [j  ]  *  1 . 1  ;minOut[j  ]  =Output  [i]  [j]  *0. 9 ; } 
else{ 

if(maxOut  [j  ]  <Output  [i]  [j  ]  *  1 . 1 )  {  maxOut  [j  ]  =Output  [i  ]  [j  ]  *  1 . 1 ; } 

if(minOut[j  ]>Output[i]  [j  ]  *  0. 9)  {  minOut[j  ]  =Output[i]  [j]  *  0. 9 ; } 

} 

} 


} 

for(i=0;i<N;i++){ 

for(j=0;j<I;j++){SealedInput[i][j]=Input_temp[i]|j]=(Input[i][j]- 
minln[j])/(maxln|j]-minln|j]); } 

for(j =0  ;j  <0;j  ++)  {  S  caledOutput  [i]  [j  ]  =Output_temp  [i]  [j  ] =(Output  [i]  [j  ]  - 
minOut  [j  ]  )/(maxOut  [j  ]  -minOut  [j  ] ); } 

} 

fcloseall(); 

} 

int  ResetWeights(){ 
int  i,j; 
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if(Run==0  ||  Run==3){//3->train  with  given  H  from  zero 


for(i=0;i<I+l;i++){for(j=0;j<H;j++){w_temp[i][j]=float(rand()%100)/100- 
0.5;w_test[i][j]=l;}}  //Small  random  initial  weights 

for(i=0  ;i<H+ 1  ;i++)  {for(j  =0;j  <0;j++)  {  vtemp  [i]  [j  ]=float(rand()%  1 00)/ 1 00- 
0.5;v_test[i]  [j]=l ;}  } 

} 

else{//Load  old  weights  from  wo.dat  &  vo.dat 
float  read; 

FILE  *fp=fopen(wo_dat,"r+"); 
fscanf(fp,"%f',&read); 

if(read!=I){printf("Error  in  weight  file  compatibility  retrain  from 
0.\n");retum  0;} 

fscanf(fp,"%i",&H); 

for(i=0;i<I+l;i++){for(j=0;j<H;j++){fseanf(fp,"%f',&w_temp[i][j]); 

w_test[i][j]=l;}} 

if(Run==2){//Read  in  Min/Max  values  used  in  training. 
for(i=0;i<I;i++){fscanf(fp,"%f',&maxIn[i]);} 
for(i=0;i<I;i++){fscanf(fp,"%f',&minIn[i]);} 

} 

fclose(fp); 

fp=fopen(vo_dat,  "r+"  ); 
fscanf(fp,"%f',&read); 

if(read!=0){printf("Error  in  weight  file  compatibility  retrain  from 
0.\n");retum  0;} 

fscanf(fp,"%i",&H); 

H_temp=H; 

for(i=0;i<H+l;i++){for(j=0;j<O;j++){fscanf(fp,"%f',&v_temp[i][j]); 
vtest  [i]  [j  ]= 1 ; }  } 

if(Run==2){//Read  in  Min/Max  values  used  in  training. 
for(i=0;i<O;i++){fscanf(fp,"%f',&maxOut[i]);} 
for(i=0;i<O;i++){fscanf(fp,"%f',&minOut[i]);} 

} 

fclose(fp); 

for(i=0;i<N;i++){ 

for(j=0;j<I;j++){SealedInput[i][j]=Input_temp[i]|j]=(Input[i][j]- 
minln  [j  ] )/ (maxln  [j  ]  -minln  [j  ] ); } 

for(j=0;j<O;j++){ScaledOutput[i][j]=Output_temp[i][j]=(Output[i][j]- 
minOut  [j  ]  )/(maxOut  [j  ]  -minOut  [j  ] ); } 

} 

} 

return  1; 

} 
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float  g(float  x,  int  i)  {//Hidden  Layer 

if(Activation==0) {return  exp(-x*x);} 
return  x*x; 

} 

float  gprime(float  x,  int  i){ 

if(Activation==0)  {return  exp(-x*  x)*  -2*  x; } 
return  2*x; 

} 

float  f(float  x){//Output  Layer 

if(Activation==0) {return  exp(-x*x);} 
return  x; 

} 

float  fprime(float  x){ 

if(Activation==0)  {return  exp(-x*  x)*  -2*  x; } 
return  1; 

} 

float  zio(int  i,  int  set){//Input  to  the  ith  hidden  neuron  assuming  all  temps  using  the 
(set)th  training  values 
float  z=0; 
int  m; 

z=w_temp[0][i]; 

for(m=0;m<I_temp;m++){z+=w_temp[m+ 1  ]  [i]  *  Inputtemp  [set]  [m] ; } 
return  z; 

} 

float  yio(int  i,  int  set)  {//Input  to  the  ith  output  neuron  assuming  all  temps  using  the  (set)th 
training  values 

float  y=0; 
int  m; 

y=v_temp[0][i]; 

for(m=0  ;m<H_temp;m++)  {  y+=v_temp  [m+ 1]  [i]  *  g(zio(m,set),m); } 
return  y; 

} 

float  SE(){//Square  scaled  error  assuming  all  temps 
float  se,yi_o; 
int  o,set; 

se=0; 

for(o=0;o<O_temp;o++){ 

for(set=0;set<N_temp;set++){ 

yi_°=yio(o,set); 

se=se+pow(Output_temp[set][o]-f(yi_o),2.0); 

} 
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} 


} 

return  0.5*se; 


void  train(){//Standard  training  with  weight  removal  (_test[]  []) 
int  maxTrainDuration=5000,  t,i,j,o,h,set; 
float  se,seo,seb,der,dertemp,dero; 
float  yi_o; 

float  norm,x,xo,dx,ds; 

float  dw  [maxl+ 1]  [maxH]  ,dv  [maxH+ 1  ]  [maxO] ; 

seo=se=seb=SE(); 

for(t=0;t<maxTrainDuration;t++){ 

norm=0; 

for(o=0;o<O_temp;o++){ 

for(h=0;h<=H_temp  ;h++)  { 
der=0; 

if(v_test[h]  [o]==l){ 
dero=l; 

for(set=0;set<N_temp;set++){ 

if(h!  =0)  {  dero=f(zio(h- 1  ,set)); } 
yi_°=yi°(°>set); 
der=der-(Output_temp  [set]  [o]- 

f(yi_o))*fprime(yi_o)  *dero ; 

} 

} 

dv[h][o]=-der;//Perform  update  all  at  once  to  avoid  error 

growth 

norm=norm+pow(der,  2. 0); 

} 

} 

for(h=0  ;h<H_temp  ;h++)  { 

for(i=0;i<=l_temp;i++){ 

der=0; 

if(w_test[i]  [h]==l){ 
dero=l; 

for(set=0;set<N_temp;set++){ 

if(i !  =0)  {  dero=Input_temp  [set]  [i- 1] ; } 
dertemp=0; 

for(o=0;o<O_temp;o++){ 

yi_°=yi°(°>set); 

dertemp=dertemp- 

(Output_temp[set]  [o]-f(yi_o))*fprime(yi_o)*v_temp[h+l]  [o] ; 

} 

der=der+dertemp*gprime(zio(h,set),h)*dero; 
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} 


} 

dw[i][h]=-der; 
norm=norm+pow(der,  2. 0); 

} 

} 

ds=pow(norm,0. 5); 

norm=l; 

x=0; 

xo=0; 

dx=mu; 

for(j =0  ;j  <  1 000 ;j ++)  { 

for(o=0;o<O_temp;o++){for(h=0;h<=H_temp;h++){v_temp[h][o]=v_temp[h]  [o] 
+dx*  dv[h]  [o]/norm; }  } 

for(h=0;h<H_temp;h++){for(i=0;i<=I_temp;i++){w_temp[i][h]=w_temp[i][h]+dx 

*  dw[i]  [h]/norm; }  } 

x=x+dx; 

se=SE(); 

if(se<=seo){seo=se;xo=x;} 

if(se>seo){ 

dx=(xo-x)/2.0; 

} 

if(fabs(dx)<(mu/128)){break;} 

} 

printf("%i\tsef=%e\txf=%f\n",t,seo,xo); 

if(se>seo){ 

dx=xo-x; 

for(o=0;o<O_temp;o++){for(h=0;h<=H_temp;h++){v_temp[h][o]=v_temp[h]  [o] 
+dx*  dv[h]  [o]/norm; }  } 

for(h=0;h<H_ternp;h++){for(i=0;i<=I_temp;i++){w_temp[i][h]=w_temp[i][h]+dx 

*  dw[i]  [h]/norm; }  } 

x=x+dx; 

} 

if(fabs(xo)<=0  &&  t>=50) {break;} 

} 

} 

float  PRtrain(){//V-fold  cross  validation  training 
float  CV,CVj; 

int  v_div=3,  //Larger  v  div  will  train  slower  but  form  better  PR 
i,j,h,o,Nd,subgroup,  NtempPR; 
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train();  //Train  with  all  Ntemp  data  sets 
CV=0; 

N_tempPR=N_temp; 

for(o=0;o<O_temp;o++){for(h=0;h<=H_temp;h++)  {v[h]  [o]=v_temp  [h]  [o] ; }  } 
for(h=0;h<H_temp;h++){for(i=0;i<=I_temp;i++){w[i]  [h]=w_temp  [i]  [h] ; }  } 

N  d=N_temp/v_di  v; 
if(v_div*Nd<N_temp){Nd=Nd+ 1 ; } 
for(subgroup=0;subgroup<N  d;subgroup++)  { 

for(o=0;o<O_temp;o++){for(h=0;h<=H_temp;h++){v_temp[h][o]=v[h][o];}}//St 
art  from  same  location  as  full  training 

for(h=0;h<H_temp;h++){for(i=0;i<=I_temp;i++){w_temp[i]  [h]=w[i]  [h]; }  } 
N_temp=0; 

for(i=0;i<N_tempPR;i++){ 

if(i<subgroup*Nd  ||  i>=(subgroup+l)*Nd){ 

for(j=0;j<O_temp;j++){Output_temp[N_temp][j]=ScaledOutput[i][j];} 

for(j=0  ;j<I_temp;j++)  { Inputtemp  [Ntemp]  [j] =ScaledInput[i]  [j  ] ; } 
N_temp=N_temp+ 1 ; 

} 

} 

train(); 

CVj=0; 

for(i=0;i<N_tempPR;i++){ 

if(i>=subgroup*Nd  &&  i<(subgroup+l)*Nd){ 

for(j=0  ;j<Otemp  ;j  ++)  {  Outputtemp  [i]  [j  ]=ScaledOutput[i]  [j] ; } 

for(j=0  ;j<I_temp;j++)  { Inputtemp  [i]  [j] =S  caledlnput[i]  [j  ] ; } 
for(j =0  ;j  <Otemp;j  ++)  {  C  Vj  =C  Vj  +po  w(Output_temp  [i][j]- 

f(yio(j,  i)),2. 0)/N_temp ; } 

} 

} 

CV=CV+CVj/v_div; 

} 

N_temp=N_tempPR; 

for(i=0;i<N_temp;i++){ 

for(j=0  ;j<Otemp  ;j  ++)  {  Outputtemp  [i]  [j  ]=ScaledOutput[i]  [j] ; } 

for(j=0  ;j<I_temp;j++)  { Inputtemp  [i]  [j] =S  caledlnput[i]  [j  ] ; } 

} 

for(o=0;o<O_temp;o++){for(h=0;h<=H_temp;h++){v_temp[h]  [o]=v[h]  [o]; }  } 
for(h=0;h<H_temp;h++){for(i=0;i<=I_temp;i++){w_temp[i]  [h]=w[i][h]; }  } 
return  CV; 
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int  SNC(){//Sequential  network  construction 
FILE  *fp=fopen("H-PR.dat","w+"); 
int  optimalH; 
int  i,j; 
float  PR; 
float  minPR; 

printf("SNC...\n"); 

I_temp=I;  0_temp=0;  N_temp=N,  H_temp=3; 
if(ResetWeights()==0){retum  -1;} 

PR=minPR=PRtrain(); 
optimalH=H_temp ; 

printf("PR(H=%i)=%f\n",H_temp,PR); 

fprintf(fp,"PR(H=%i)=%f\n",H_temp,PR); 

for(H_temp=5;H_temp<=H;H_temp+=2){ 

for(i=0;i<I  temp+l;i++){for(i=0;i<H  temp-l;i++){w  testriiril=0;}} 
for(i=0;i<H_temp;i++){for(j=0;j<O;j++){v_test[i][j]=0;}} 
train();//lst  train  with  old  weights  held  fixed 

for(i=0;i<I  temp+l;i++){for(i=0;i<H  temp-l;i++){w  test[i][il=l;}} 
for(i=0;i<H_temp;i++){for(j=0;j<O;j++){v_test[i][j]=l;}} 
PR=PRtrain();//2nd  train  on  all  weights 
if(PR<minPR){ 

minPR=PR;  optimalH=H_temp; 

for(i=0;i<I_temp+l;i++){for(j=0;j<H_temp;j++){w_optimal[i][j]=w_temp[i][j];} 


for(i=0  ;i<=H_temp;i++){for(j  =0;j  <Otemp  ;j++)  {  voptimal  [i]  [j  ]  =v_temp  [i]  [j  ] ; }  } 

} 

printf("PR(H=%i)=%f\tSE=%f\n",H_temp,PR,SE()); 

fprintf(fp,"PR(H=%i)=%f\tSE=%f\n",H_temp,PR,SE()); 

if(PR<0.0000005){break;} 

} 

for(i=0;i<I_temp+l;i++){for(j=0;j<optimalH;j++){w_temp[i][j]=w_optimal[i][j]; 

}} 

for(i=0;i<=optimalH;i++){for(j=0;j<O_temp;j++){v_temp[i]|j]=v_optimal[i][j];} 

} 

fcloseall(); 
return  optimalH; 

} 

void  PrintFit(){ 

FILE  *fp=fopen("Fit.dat",Mw+"); 
int  o,  set; 
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float  yn,yr; 

fprintf(fp,"  H=%i\n"  ,H_temp); 
fprintf(fp,"True\tFit\tError\n"); 
printf("H=%i\n",H_temp); 
printf("True\tF  it\tError\n"  ); 
for(set=0;set<N_temp;set++){ 

for(o=0;o<O_temp;o++){ 

yn=f(yio(o,set))*(maxOut[o]-minOut[o])+minOut[o]; 
yr=ScaledOutput[set][o]*(maxOut[o]-minOut[o])+minOut[o]; 
fprintf(fp,"%f\t%f\t%f\n",yr,yn,fabs(yr-yn)/(yr)*  100); 
printf("%f\t%f\t%f\n",yr,yn,fabs(yr-yn)/(yr)*  1 00); 

} 

} 

fcloseall(); 

} 

void  PrintWeights(){ 
int  i,j; 

FILE  *fp=fopen(wo_dat,"w+"); 
fprintf(fp,"%i\n%i\n",I_temp,H_temp); 

for(i=0;i<I_temp+l;i++){for(j=0;j<H_temp;j++){fprintf(fp,"%f\t",w_temp[i][j]);} 

fprintf(fp,"\n");} 

fprintf(fp,"\n"); 

for(i=0  ;i<I_temp;  i++)  {fprintf(fp,"  %f\t",maxln  [i]) ; } 
fprintf(fp,"\n"); 

for(i=0  ;i<I_temp;  i++)  {fprintf(fp,"  %f\t",minln[i] ); } 
fclose(fp); 

fp=fopen(vo_dat,  "w+"); 

fprintf(fp,"  %i\n%i\n"  ,Otemp,  Htemp); 

for(i=0  ;i<H_temp+ 1  ;i++)  {for(j  =0  ;j<0;j++)  {fprintf(fp,"  %f\t",v_temp  [i]  [j  ] );  }fprint 
f(fp,"\n");} 

fprintf(fp,"\n"); 

for(i=0;i<O;i++){fprintf(fp,"%f\t",maxOut[i]);} 

fprintf(fp,"\n"); 

for(i=0;i<O;i++){fprintf(fp,"%f\t",minOut[i]);} 

fclose(fp); 


int  main(){ 

LoadData(); 

if(Run==2){ 

ResetWeights(); 

H_temp=H;  I_temp=I;  0_temp=0;  N_temp=N; 
PrintFit(); 

} 
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else{ 


if(Run==0){ 

H_temp=SNC(); 

I_temp=I;  0_temp=0;  N_temp=N; 
train(); 

if(H_temp==-l){retum  0;} 

} 

else{ 

ResetWeights(); 

H_temp=H;  I_temp=I;  0_temp=0;  N_temp=N; 
printf("Re-training  with  H=%i  from  last  training 

values\n",H_temp); 

train(); 

} 

PrintWeights(); 

PrintFit(); 

} 

return  1; 

} 
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